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1 . INTRODUCTION 


This set of notes discusses the numerical solution of transonic flow 
fields governed by the full-potential equation. In a general sense this 
presentation deals with relaxation schemes suitable for the numerical solu- 
tion of elliptic partial differential equations. Of course, transonic flow 
is not purely elliptic in nature but consists of hyperbolic regions embedded 
In an otherwise elliptic domain. However, the most successful numerical 
methods of solution for transonic flow applications, at least for potential 
formulations, have evolved from classical relaxation schemes associated with 
elliptic equations. Thus, most of the algorithms presented herein will have 
an elliptic-equation, relaxation-algorithm flavor. For related material the 
reader is referred to Hall (ref. 1) for a review with emphasis on the his- 
torical aspects of computational fluid dynamics (CFD) , Holst et al. (ref. 2) 
for a general review of computational transonic aerodynamics (CTA) , Chapman 
(refs. 3,4) for a general review of current CFD research, and Kordulla 
(ref. 5) and Baker (ref. 6) for additional information regarding transonic- 
flow solution methods. 

The transonic speed regime provides the most efficient aircraft cruise 
performance; hence, most large commercial aircraft cruise in this speed 
regime. However, transonic flow fields tend to be sensitive to small per- 
turbations in flow conditions or to slight changes in geometrical character- 
istics, either of which can cause large variations in the flow field. Large 
performance penalties can result because of relatively small perturbations 
away from desired design conditions. Computational techniques, therefore, 
have enjoyed an increasing role in helping the aerodynamics engineer find 
optimal design conditions, as well as to evaluate design sensitivity. 

Transonic flow fields contain a variety of interesting and unique char- 
acteristics. Typical airfoil and swept-wing flow fields are shown in 
figures 1 and 2. The outer free-stream flow is typically subsonic and 
elliptic in nature. Regions of supersonic flow (hyperbolic) usually exist 
on the upper airfoil or wing surface and are generally terminated by a weak 
’'transonic" shock wave. For the case of a swept-wing flow field, the shock 
wave may actually consist of a system of shocks, as shown in figure 2. The 
first shock is swept and therefore has a supersonic downstream Mach number. 
The rear shock is approximately normal to the local flow and therefore has 
a subsonic downstream Mach number. 

Signals tend to propagate very rapidly downstream in transonic flow 
fields where the propagation velocity is u + a (local fluid speed plus sonic 
velocity) and very slowly upstream where the propagation velocity is u - a. 
For a downstream disturbance to propagate upstream it must move around the 
supersonic zone, further increasing the difference between the upstream and 
downstream propagation speeds. This situation tends to make transonic 
numerical solution techniques, which depend on physical-time-dependent algo- 
rithms, very slow. Such problems are said to be "stiff," and they require 
large amounts of computer time for even a single steady-state solution. 
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Another characteristic of transonic flow is that it is governed by 
equations that are inherently nonlinear. Linearization of these equations 
will remove the vital flow-field physics which are responsible for the pre- 
diction of shock waves. This nonlinear behavior tells us that a direct 
(i.e., noniterative) solution procedure for transonic flow is theoretically 
impossible. Thus, one basic feature associated with all of the schemes dis- 
cussed within these notes is that they are all iterative. 

Viscous effects are also important in transonic flow fields. This 
complex subject involves four major effects: (1) shock/boundary-layer inter- 

action effects, (2) the decambering and thickness effects caused by the 
addition of a simple displacement thickness, (3) trailing-edge effects, and 
(4) near-wake effects. However, a discussion of viscous correction proce- 
dures is not within the scope of the present set of notes. Instead, the 
interested reader is referred to Lock (ref. 7) and Melnik (ref. 8) who treat 
this subject in detail. 

The following presentation can be divided into two general categories: 
an introductory treatment of the basic concepts associated with the numerical 
solution of elliptic partial differential equations and a more advanced 
treatment of current procedures used to solve the full-potential equation 
for transonic flow fields. The introductory material is presented for com- 
pleteness; it covers governing equations (sec. 2), classical relaxation 
schemes (sec. 3), and early concepts regarding transonic, full-potential 
equation algorithms (sec. 4). These topics are intended to provide an intro- 
duction for some of the more advanced concepts presented in the later sec- 
tions. The more knowledgeable reader could skip sections 2-4 and proceed 
directly to sections 5-8 without a significant loss in continuity. 

State-of-the-art topics concerning the numerical solution of the full- 
potential equation for transonic flows are presented in sections 5-8. These 
sections include a discussion of equation transformation and grid-generation 
procedures (sec. 5); a presentation of recent full-potential spatial differ- 
encing schemes (sec. 6); a presentation of full-potential iteration schemes, 
with special emphasis on convergence acceleration (sec. 7); and a brief 
review of recent three-dimensional applications (sec. 8). These notes are 
then concluded with a few general remarks including recommendations for 
future research. 


2. GOVERNING EQUATIONS 


2.1 Classification of Second-Order Partial Differential Equations 

Consider the following general quasi-linear , second-order partial dif- 
ferential equation (PDE): 


Au 


XX 


Bu + Cu 
xy yy 


= F 


( 2 . 1 ) 


where A, B, C, and F are functions of x, y, u, u^, Uy. This equation can 
be classified by considering the corresponding characteristic equation (for a 
derivation of the characteristic equation and additional discussion on this 
topic see refs. 9 and 10): 



- B ^ + C = 0 
dx 


( 2 . 2 ) 
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Using the quadratic formula, the two characteristic directions associated 
with equation (2.1) are given by 


dy\ B ± /B^ - 4Ac 
‘^"4,2 “ 2A 


(2.3) 


The nature of these characteristics determines the equation classification. 
Equation (2.1) is hyperbolic If the characteristics are real and distinct 
[that is, If the discriminant of eq. (2.3) Is greater than zero 
(B^ - 4AC > 0), then the equation is hyperbolic]; parabolic if the character- 
istics are real and coincidental (B^ - 4AC = 0); and elliptic if the charac- 
teristics are complex and distinct (B^ - 4AC < 0). 

Several classical FDEs along with their classifications are given as 
follows . 

Laplace’s equation: 


u 


XX 


+ u 


yy 


0 


Heat equation: 


Wave equation: 


= - 

elliptic 

u = a^u 
y XX 

(o ~ real) 

= - 

parabolic 

u = c^u 

(c ~ real) 

XX yy 


(n. = - 

hyperbolic 


The primary motivation for studying the nature of PDEs in the present 
context is to gain insight into the physics they govern, and, therefore, to 
develop guidelines for the implementation of numerical solution procedures. 
Different equation types generally require different solution algorithms. 


A particular point, P, in a solution domain governed by a PDE, has 
associated with it regions called domains of influence and domains of depen- 
dence. These domains are determined by characteristics. For example, in 
the case of steady supersonic flow, which is hyperbolic, the domain of depen- 
dence is defined by the forward Mach cone, and the domain of influence is 
defined by the aft Mach cone (see fig. 3). Physically, the flow conditions 
at point P depend on the flow conditions in the domain of dependence. Con- 
versely, the flow conditions at point P can influence any other point in 
the domain of influence. 

The numerical domain of dependence, that is, the domain of dependence 
modeled by the numerical differencing scheme, should mimic the physical 
domain of dependence as closely as possible. If it does not, instability 
probably will result. 
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Elliptic equations are much different in nature than hyperbolic equa- 
tions. A single point in an elliptic solution domain influences every other 
point and vice versa; that is, the elliptic domains of dependence and influ- 
ence are coincidental and contain the entire solution domain. Solution 
algorithms for elliptic equations should reflect this fact to properly 
simulate the physics. 

For transonic flow applications, the solution domain contains both 
elliptic (subsonic) and hyperbolic (supersonic) regions. The boundaries 
between these regions (sonic lines or shock waves) are unknown in advance 
and must be determined as part of the solution. Most successful relaxation 
schemes used for transonic applications utilize a feature called type- 
dependent differencing which allows the local differencing scheme to be 
adjusted to the local flow-field physics. Numerical solution techniques for 
implementing this philosophy will be the primary subject of discussion in 
subsequent sections of these notes. Two PDE formulations representing 
different levels of approximation for transonic flow fields are discussed 
next. 


2.2 Transonic, Small-Disturbance Potential Equation 

The transonic, small-disturbance (TSD) equation expressed in two- 
dimensional Cartesian coordinates (x,y) Is given by 

[1 - - m2(y + + <^yy = 0 (2.4) . 

where is the free-stream Mach number, y is the ratio of specific heats, 

and <f> is the small-disturbance or perturbation velocity potential defined 
by 

V<P (2.5) 

In equation (2.5), q and q^ are the local and free-stream velocity vectors 
defined by 

q = ui -f vj , q = Uj. (2.6) 

where i and j are the unit vectors along the x and y directions, respec- 
tively. Thus, the velocity components are given by 

= U - , (J> = v (2.7) 

X y 

Boundary conditions for a typical "thin" airfoil used in conjunction 
with equation (2.4) are given by 

+ 

4>y(x,0-) = U„ (2.8) 

where g^(x) and g”(x) define the upper and lower airfoil surfaces, respec- 
tively. Application of this boundary condition is made at the airfoil slit 
(or chord line) . This simulates the required f low-tangency boundary condi- 
tion at the airfoil surface to an accuracy consistent with small-disturbance 
theory. An auxiliary relation, usually used in this formulation to define 
the airfoil surface pressure coefficient, is given by 
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cj = -2(|.^(x.O-) (2.9) 

where Cp and Cp correspond to the upper and lower surfaces, respectively. 

The TSD equation, given by equation (2.4), is valid for Isentroplc, 
Irrotational flows about thin shapes (airfoils and wings) Immersed in approx- 
imately sonic free streams (that is, ~ 1). For a more detailed discussion 
of TSD theory, including a discussion of the various three-dimensional formu- 
lations, see references 2 and 11. 

By using the discriminant test obtained from equation (2.3), it can be 
shown that equation (2.4) is hyperbolic when 

^ . 

(y + ml 

and elliptic when 

(j, < T Hw — 

(y + ml 

In other words, the sign of the first term coefficient determines the equa- 
tion type. If this coefficient is positive, the local flow is subsonic; if 
it is negative, the local flow is supersonic. The nonlinearity of the first 
term is essential for describing the mixed character of transonic flow and 
is the mechanism by which shock waves are formed. 

The characteristic directions associated with the TSD equation are given 
by 

(S) = ±[1 - Mi - Mi(Y + 


Notice that these characteristic slopes are symmetrical about the x-axis 
regardless of the local velocity-vector orientation; that is, the character- 
istics are not a function of the y-component of the velocity (i|)y) (see 
fig. 4). This situation, which is in dramatic contrast to the full-potential 
equation, has certain implications regarding spatial-difference approxima- 
tions for both the TSD and the full-potential equation. 


2.3 Full-Potential Equation 


The full-potential equation written in conservative form for two- 
dimensional Cartesian coordinates (x,y) is given by 

(P<t>^)x + (P'^y)y = 0 (2.11a) 

where the density p is defined by 

In equation (2.11), y is the ratio of specific heats (equal to 1.4 for air) 
and (|) is the full or exact velocity potential given by 

Vf|i = <fi i + 4> j = ui + vj = q (2.12) 

X y 
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where ijix and ifiy are the velocity components in the x and y directions, 
respectively. The density and velocity components appearing in equa- 
tion (2.11) are nondimenslonalized by the stagnation density Pg and by the 
critical speed of sound, a^. Thus, at a stagnation point. 

‘*’x*'^y ® * P 1 

(2.13) 

and at a sonic line. 


<t>J + 1 . P [l - “ 0.633938145 

. . . (2.14) 

The latter condition is quite useful in providing a simple test for super- 
sonic flow; that is, if p < 0.633938145 . . . the flow is supersonic, and 
if p > 0.633938145 . . . the flow is subsonic. 

Several useful auxiliary relations are given as follows. 


Bernoulli equation (energy equation): 


2 2 
Q ^ a ^ 

+ r = constant 

2 Y - 1 

(2.15) 

Isentropic equation of state (perfect gas): 


P/p^ = constant 

(2.16) 

Speed of sound definition: 


.2 = dP _ IP 

dp p 

(2.17) 

The second equality in equation (2.17) is obtained by using the isentropic 
equation of state. Using the same nondimenslonalizatlon for these realtlons 
as for the full-potential equation yields 

a^ 1 Y + 1 

2 Y - 1 2 Y - 1 

(2.18) 

P Y + 1 

pY " 2 y 

(2.19) 

a^ = li 

n 

(2.20) 


Note that the density expression shown in equation (2 .lib) has been derived 
from equations (2 . 18)-(2. 20) . 

The full-potential equation given by equation (2.11) is valid for isen- 
tropic, irrotational flows about arbitrary shapes. That is, the full- 
potential equation is not restricted to thin shapes as is the TSD equation. 
However, to obtain physically realistic results, the full-potential equation 
is restricted to shapes and to flows for which viscous effects (separated 
flow) are not significant. The full-potential equation is also restricted 
to flows ranging from incompressible (M<„ ~ 0) to transonic (M<o ~ 1) • These 
potential formulations are valid in transonic flow (even though they are, 
strictly speaking, isentropic) because weak, or transonic, shock waves can be 
reasonably approximated by an isentropic formulation. A comparison of the 
isentropic shock-jump relation with the exact Euler shock-jump relations 
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(Rankine-Hugoniot relations) Is given in figure 5 (see ref. 12), Note that 
for a local Mach number (U-^) at or below 1.3 a reasonable approximation is 
maintained by the Isehtroplc assumption. 

It is essential that the finite-difference approximation to the full- 
potential equation be cast in conservative form (see refs. 13 and 14). 
Otherwl.se • the shock-capturing procedure will not necessarily conserve mass 
across the shock wave. The effective mass source added at the shock wave Is 
dependent on nonphysical considerations such as the local grid spacing. 
Nonconservative schemes have been used, with good success, in transonic-flow 
simulations for many engineering applications. This is because the effec- 
tive mass production at shocks fortuitously models shock/boundary-layer 
interactions. A superior approach is to use conservative form with viscous 
corrections added via the boundary- layer equations. Therefore, the proper 
shock strength and position (within the framework of the isentroplc formula- 
tion), automatically corrected by the proper viscous effects, will be 
obtained. 

Classification of the full-potential equation [eq.. (2.11)] Is difficult 
because it is not in standard form [eq. (2.1)]. Transformation of equa- 
tion (2.11) into nonconservative form will facilitate the classification 
process. The nonconservative full-potential equation written in two- 
dimensional Cartesian coordinates is given by 

where a is the local speed of sound computed from the Bernoulli equation 
[see eq. (2.18)]. The discriminant [see eq. (2.3)] for equation (2.21) 
becomes 


- 4AC =« 4a^(q^ - a^) (2.22) 

Therefore, the full-potential equation, either in the conservative form 
[eq. (2.11)] or nonconservative form [eq. (2.21)], is hyperbolic for super- 
sonic flow (q > a) , parabolic for sonic flow (q = a) , and elliptic for sub- 
sonic flow (q < a). 


The characteristic directions associated with the full-potential equa- 
tion are given by 


( 


,dx/^ 


a" - 


(2.23) 


Notice that the characteristic directions for the full-potential equation are 
not symmetric about the x-axis as was the case with the TSD equation. 
Instead, the characteristics are symmetric about the stream direction. (This 
fact will become obvious in a subsequent section.) Therefore, like the local 
velocity vector q, the characteristic directions can have virtually any 
orientation, which makes numerical solution of the full-potential equation 
more complicated than the previously discussed TSD equation. After the pres- 
entation of several classical relaxation schemes for purely elliptic equa- 
tions, we shall return to this and other points associated with the numerical 
solution of the full-potential equation for transonic flows. 
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3. CLASSICAL RELAXATION ALGORITHMS 


A look at some of the classical relaxation schemes applied to elliptic 
equations, including point and block (line) iterative techniques, Is consid- 
ered next. Additional details on this subject can be found in Ames (ref. 9). 
To facilitate this discussion, notation conventions, which will apply 
throughout these notes, are now presented. 


3.1 Notation Conventions 

In general, the notation, will be used to represent the nth 

iterate of the discretized dependent variable, ij), at a position in the 
finite-difference mesh given by x = lAx and y * jAy. When transformed 
coordinates (?,n) are involved, the i and j subscripts will be used to 
represent position in the computational-domain, finite-difference mesh (that 
■ is, ? “ i and n = j where A? = An = 1.0). Definitions for the various 
difference operators used in these notes are given by the following. 

First central difference (second-order accurate): 


^x^ ^i.j = 2Ax ^i+i,j “ ^ ^i-i,j^ 
First backward difference (first-order accurate): 

^x^ ^i.j = 'h. ^i,j ■ ^ ^i-i,j’ 
First forward difference (first-order accurate): 

^i,j - Ax ^i+i,j " ^ ^i,j^ 
Second central difference (second-order accurate): 

> 1.3 + < 

Forward shift operator: 


Backward shift operator: 


Central average operator: 


^x ^ ^i.j = ^ ^1+1. j 


>i.j = ( >i-x,j 


_ 1 


^x^ ^i,j " 2 ^i+i,j ^ ^i-i,j^ 

Forward average operator: 


^i,j " 2 ^i+i,j ^ ^i.j^ 


(3.1) 

(3.2) 

(3.3) 
(3.A) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 
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Backward average operator: 






(3.9) 


These operator definitions have all been defined using the x~dlrectlon. Of 
course, additional operators can be defined for any space or time coordinate 
similarly to those for the x-dlrection. Certain Identities associated with 
these operators exist and can be useful. For example, all difference and 
averaging operators can be expressed, using shift-operator notation: 




n, - i (I + 


t t - -L (1 _ e"^) (1 - e“^) 

XX Ax X Ax X 


E"^6 

X XX 


(3.10) 


= — (1 - 2E ^ + E ^) 

A 2 * ^ 

Ax 

Notice that the symmetrical combination of first-order-accurate, one-sided 
difference operators creates second-order-accurate, centered-dlfference 
operators. 


1 +1 -1 = x ' - ^ 1 

6 = -— (E^ - 2 + E ^) - r-^ 

XX , 2 X X Ax Ax 

Ax 


XX XX 


6 = ~ (e"^^ - E"^) » (■? + t ) 

X 2Ax X x' 2 x x' 


(3.11) 


This fact will be useful in the application of fully Implicit, approximate- 
factorization schemes for solving the full-potential equation. 

The classical relaxation schemes discussed next will all be presented 
for the two-dimensional Laplace equation: 


XX 


(b 

‘yy 


(3.12) 


Each scheme will be put into a standard, two-level correction form (sometimes 
called delta form) given by 


NC" , + (uL(}>" . = 0 


(3.13) 


where 


c" .. is the 
1.3 


nth-iteration correction defined by 


,n _ ,n+i 
'l.j ■ 


(3.14) 


L((i2 a is the nth-iteration residual, which is a measure of how well the 
finite-difference equation Is satisfied by the nth-level solution, and 
0 ) is a relaxation parameter. The general Iteration scheme given by 
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equation (3.13) can be considered as an Iteration in time (that iSt pseudo- 
time — the scheme does not actually apply to a physical time-dependent dif- 
ferential equation) . This consideration allows the n superscript to be 
regarded aa a time-index, that is, ~ At(t>t« N-operator determines the 

type of iteration scheme and, therefore, is the or.ly quantity from equa- 
tion (3.13) to change from scheme to scheme. 


3.2 The von Neumann Tbst for Stability 

The general iteration scheme represented by equation (3.13) can be 
investigated for numerical stability by using the Fourier, or von Neumann, 
teat for numerical stability c This scheme was developed by von Neumann in 
the early 1940s and was first discussed in detail by O'Brien et al. (ref. 15) 
in 1951. Additional details can be found in Smith (ref. 16), Ames (ref. 9), 
Mitchell (ref. 10), and Richtmyer and Morton (ref. 17). 

The propagation of numerical error is studied by substituting a suitable 
solution. 



Xt iax iby 
e e e •' 


(3.15) 


into the finite-difference scheme [eq. (3.13)]. In equation (3.15) a and b 
are wave numbers associated with the Fourier series which represents the 
solution at t «> 0, X is (in general) a complex constant, and the super- 
script i is /^. The resulting expression is solved for 

G = g>^(t+At)^gXt ^ gXAt 


This quantity is called the amplification factor and provides an indication 
of error growth or decay. If G > 1 for any values of a and b, the scheme 
is numerically unstable. If G < 1 for all values of a and b, the scheme 
is stable. 


There are several important limitations associated with the von Neumann 
stability test. First, this test applies only to linear difference schemes 
with constant coefficients; however, if the difference scheme in question has 
variable coefficients, the method can still be applied locally with a good 
chance for accurately predicting numerical stability characteristics. 

Second, the effect of boundary conditions is neglected by this method (that 
is, the boundary conditions are assumed to be periodic). Another type of 
stability analysis, referred to as the matrix method, is preferable when the 
effect of boundary conditions needs to be analyzed (see refs. 9, 10, 16 for 
discussion on this topic) . 


3.3 The Point-Jacobl Scheme 

The point-Jacobi relaxation scheme can be expressed as follows: 


/a" 




Rewriting this scheme in standard form yields 


/ _ 2 _ . _2_Vn 

\ Ax^ Ay 7 


+ )i|). , 

XX yy'^lsj 


(3.18) 


10 


where 


N • N 


PJ 


Ax" 


Ay-* 


and 


L 


6 +6 
XX yy 


(3.19) 


(3.20) 


Note that the N-operator for the polnt-Jacobl scheme Is Just a scalar or In 
matrix form a simple diagonal matrix. Notice also that u ■ 1; that is, the 
point-Jacobi scheme is not overrelaxated. Decause of the simple form of 
Npj, obtaining the (n + l)st-level solution from the nth-level solution is 
quite easily accomplished. For Instance, rewriting equation (3.17) with 
only the n + 1 term on the left and with Ax •• Ay yields 



^^i+i,j ^ Vi. 








(3.21) 


Thus • the value at each grid point is simply replaced by the average of the 
values at the four nearest neighboring points. The simplistic nature of the 
point-Jacobi scheme and, in particular, the large difference between the 
construction of the N and L o>perators, yields very slow convergence, 
especially as the mesh is refined. 


Stability for the point-Jacobi scheme can be investigated by using the 
von Neumann stability analysis (see sec. 3.2). The amplification factor 
(for Ax ■ Ay) is given by 

Gpj - 1 + j u -3 aAx - 1) + I (cos bAy - 1) (3.22) 


For stability, the magnitude of the amplification factor must be bounded 
by 1. Notice that Gpj is purely real. This makes the stability analysis 

somewhat easier. Since both cos aAx - 1 and cos bAy - 1 range from a 

maximum of sero to a minimum of -2, the growth factor can never be larger 

than 1 nor smaller than -1. Notice that if u had been retained, any value 

above 1 (representing over relaxation) would have caused instability. 

If equation (3.18) is rewritten as 


♦ 


n+x 

ij 



At (6 


XX 


6 

yy 



(3.23) 


where At « (2/Ax^ + 2/Ay^)“^, the point-Jacobi relaxation scheme looks very 
much like a standard heat-equation integration scheme in which At is the 
time-integration step-size. The stability condition associated with the 
time-integration of the heat equation using the present "explicit scheme" is 
given by 


At/Ax* + At/Ay* i 1/2 (3.2A) 

This condition is automatically satisfied by the point-Jacobi scheme through 
the definition of At. 
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3.4 The Polnt"Gauss-Seidel Scheme 


The point-Gauss-Seidel relaxation scheme is very similar to point-Jacobl 
with one Important difference: as we sweep through the mesh, the solution at 

some points has already been updated. Polnt-Gauss-Seidel uses this "extra" 
information by using the latest updated solution at every point possible. 

For instance, assuming we sweep through Increasing values of the mesh Indices, 
1 and j , we have 



(<!>’ 


n 


1+1, j 


2<|) 


n+i 

i.j 


,n+i 




Ay* 




n 


i.j+i 


2-^ 


n+i . .,n+i 


i.j ' '^l.J 




0 (3.25) 


Note that the 1-1, j and l,j-l terms, as well as the usual l,j term, are 
all evaluated at the new n + 1 iteration level. Rewriting equation (3.25) 
In standard form yields 


where 


(E~' - 2 ) +-^ (E"^ - 2) C" + - 0 (3.26) 

Ax* * Ay* y J 



(3.27) 


and again w ® 1. The residual operator L is defined as before by equa- 
tion (3.20). The N-operator for the point-Gauss-Seidel scheme (Npgg) pro- 
vides a more complete approximation to the L-operator than Npj. Because of 
this, the point-Gauss-Seidel scheme convergence is faster than that of point- 
Jacobi. As we shall see, however, other schemes with much faster convergence 
speed exist. 

Because point-Gauss-Seidel has several entries per row, located on or to 
the left of the diagonal, advancing from the nth-solution level to the 
(n + l)st-solution level does not Involve typical matrix Inversions, for 
example, tridiagonal and pentadiagonal. However, obtaining the (n + l)st- 
solution level is basically a recursive process. Such recursion is difficult 
(or imposRible) to program on vector computers in such a way as to fully 
utilize the improvement in computational speed offered by the vector archi- 
tecture. In this context, point-Jacobi , which is completely explicit, that 
is, nonrecursive, may actually be more efficient than point-Gauss-Seidel. 
Other relaxation schemes which are vectorizable are discussed in section 7.6. 


3.5 The Successive Over relaxation Scheme 


The successive overrelaxation (SOR) scheme can be presented by first 
writing the point-Gauss-Seidel scheme (3i25)] in the following form: 



(<^ 


i+i,j 


<] ^ ^ 


n+i 

i-i.j 


) + 


Ay^ 


(f 


l,j+i 


2*: 


n+i 


+ <!' 


n+i 

i.j-i 


) 


0 

(3.28) 


In equation (3.28) the values of n + 1 with an over bar are provisional 
values modified by the following standard relaxation formula: 


-f- 


n+i 

i,j 



+ (1 - r)(|)J^j 


(3.29) 
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where r is a relaxation factor. For values of r > 1 the scheme is said 
to be overrelaxed, and for values of r < 1 the scheme is said to be under- 
relaxed. Equations (3.28) and (3.29) can be combined and rewritten in 
standard form as 



(3.30) 


The N-operator for the SOR scheme is still of the same form as NpQg 
[eq. (3.27)], although the matrix elements have been modified by the pres- 
ence of r. Use of the overrelaxation factor greatly improves the conver- 
gence rate of the SOR scheme relative to the polnt-Gauss-Seidel scheme, as we 
shall see in section 3.8. 


Stability of the SOR scheme as applied to Laplace'^s equation can easily 
be investigated by the von Neumsmn analysis. The amplification factor for 
the SOR scheme is given by 

4 

-cos aAx - cos bAy - — + 4 - i(sin aAx + sin bAy) 

GgoR = 4 (3.31) 

cos aAx + cos bAy ~ ~ ~ i(sin aAx + sin bAy) 


Notice that the imaginary parts of both the numerator and denominator are 
identical. Therefore, an alternative, but no less restrictive, condition for 
stability is given by 


-cos aAx - cos bAy +4 

_1 < L_ < 1 

cos aAx + cos bAy 


(3.32) 


Assume r is such, that the denominator is negative. This yields, after 
simplification, the following two Inequalities: 

2/r > 1 (3.33) 


-cos aAx - cos bAy > -2 


(3.34) 


which are satisfied for all values aAx and bAy providing 0 < r < 2. The 
latter restriction on r forces the denominator to always be negative, which 
is consistent with out initial assumption. Thus, the SOR iteration scheme is 
stable by the von Neumann stability test. 


An interesting approach for finding an optimum value for the SOR relaxa- 
tion parameter was presented by Garabddlan (ref. 18; see also ref. 9). The 
approach consisted of viewing the SOR solution procedure as a time-dependent 
problem, that is, an artificial-time-dependent problem. Major' features of 
this technique can be outlined as follows: First, using Taylor series derive 

the modified equation, then assume that Ax and At are so small that the 
solution to the difference equation (or modified equation) behaves like the 
solution to the corresponding differential equation. The time-dependent PDE 
can then be solved analytically using a separation of variables technique. 

The pvoblem is to choose a relaxation factor r such that the time-dependent 
analytical solution converges as rapidly as possible to the desired steady- 
state solution. A typical result obtained with this approach using a square 
domain (0 < x < tt and 0 < y < it) is given by 


2 

K ■ ' 

opt 1 + Ax 


(3.35) 
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Thus* the optimum relaxation factor approaches a value of 2 as the finite- 
difference mesh is refined. This analysis was performed on Laplace's equa- 
tion, and, a3,though strictly speaking it is not valid for nonlinear problems, 
the same results apply qualitatively. 


3.6 Block or Line Iteration Schemes 

Thus far only point iterative schemes have been discussed. Other 
methods which involve the simultaneous evaluation of more points at the 
n + 1 Iteration level are now considered. This additional implicitness pro- 
vides for faster convergence while requiring only minlmel increases in com- 
putational work. The first three line schemes discussed are line-Jacobi, 
line-Gauss-Seidel, and successive line-overrelaxation (SLOR). These schemes, 
written in standard correction form, are given by the following. 

Line-Jacobi : 





+ L(|) 


i.j 


0 


Line-Gauss-Seidel ; 


r_i_ 

Lax^^ 



- 2 ) + 6 




+ L(j> 


l.j 


0 


SLOR; 




-Ax 


2 V^x'“ 




+ L<j)' 


i.j 


“ 0 


(3.36) 


(3.37) 


(3.38) 


In each case one additional point has been added to the N-operator differ- 
encing molecule. This allows the entire y-dlfference operator to be modeled 
in the N-operator and requires the inversion of a set of tridiagonal matrix 
equations for each iteration. Equations (3.36)-(3. 38) have all been written 
for vertical-line relaxation, that is, the tridiagonal matrix extends along 
the y-direction. An equally valid form can be constructed by using 
horizontal-line relaxation. The optimal form depends on characteristics of 
the particular application. 

The SLOR scheme [eq. (3.38)] is widely used today in industry as a 
general relaxation scheme for solving both the TSD and the full-potential 
equations. A specific SLOR algorithm used to solve the nonconservative full- 
potential equation will be discussed in section A. 


3.7 The Alternating Direction Implicit Scheme 

One technique for achieving even faster convergence than that provided 
by the SLOR scheme is to use a fully implicit scheme; that is, a scheme in 
which each point communicates with every other point during each interatlon. 
This type of scheme can be constructed using the approximate factorization 
(AF) philosophy. First, write the N-operator in direct form, that is, 

N = L, where L is the usual residual operator given by equation (3.20). 

Note that if the problem is nonlinear, then N must be a linearized approxi- 
mation to the nonlinear L-operator. The next step in the construction of a 
suitable AF scheme is to factor the N-operator by using an appropriate 
factorization, as indicated by 
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N ~ N.Nj 


(3.39) 


Ueually the AF scheme consists of two factors for two-dimensional problems 
and three factors for three-dimensional problems. An important idea behind 
the factorization is that each of the factors Nx and Nj must involve only 
simple banded matrix inversions, thus reducing the computational v;ork per 
iteration. Then, both the errors associated with the factorization and the 
linearization (for nonlinear problems) are removed from the solution simul- 
taneously and automatically by iteration. 

Of course, using the N-operator directly (before factorization) as the 
final form of the N-operator will yield a direct scheme which converges in 
just a single iteration (for linear problems). This type of scheme requires 
a special inversion algorithm (direct solver) which is much more computa- 
tionally expensive than, for instance, a set of tridiagonal matrix inver- 
sions. For nonlinear problems, which must be considered for transonic appli- 
cations, the direct-solver scheme loses its single-iteration advantage and 
must be iterated just as any other standard scheme. Therefore, the direct- 
solver scheme will not be considered further in these notes. For more 
information see references 19-23. 

The first AF scheme presented is a reformulation of the Peaceman- 
Rachford alternating direction implicit (ADI) scheme and can be expressed by 
choosing N as follows (see refs. 9 and 10 for more information on ADI 
schemes) : 


'ADI 


(ct 


6 ) (a 

XX 




(3.40) 


where a is an acceleration parameter which may be considered as the inverse 
of an artificial time-step, l/4t. More on the optimal choice of the a 
parameter will be presented subsequently. In equation (3.40) both the x and 
y directions are treated implicitly. The N-operator has been written as 
the product of two tridiagonal matrix factors which when multiplied out yield 


ADI 


-o - 


6 (S + L 

XX yy 


(3.41) 


This expression contains the original L-operator plus two error terms. The 
ADI scheme of equation (3.40) can be restated in practical terms using two 
sweeps as follows. 

Sweep 1: 

" \x>^i.j " ^^-^2) 

Sweep 2: 

(a - 6 )g" => f" (3.43) 

yy' i.j i.j 

In equations (3.42) and ^3.43) u is the relaxation parameter associated with 
the standard form and is an intermediate result stored at each mesh 

point in the finite-difference mesh. In the first sweep, the f-array is 
obtained by solving a tridiagonal matrix equation for each y « constant 
line. The correction array (Ci,j) is then obtained in the second sweep from 
the f-array by solving a tridiagonal matrix equation for each x •• constant 
line. This construction of N allows each grid point in the entire mesh to 
be influenced by every other grid point during each Iteration. As a result. 
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much faster convergence can be obtained with this type of algorithm. More 
discussion on the ADI scheme in particular and AF schemes in general will be 
presented in section 7 . 


3.8 Convergence Rate Estimation 

A summary of the schemes presented in the preceding sections of this 
chapter is given in figure 6. The computational molecules of both the N and 
L operators are shown schematically for each scheme. Note that the 
L-operator consists of the same five grid points in each case whereas the 
N-operator is different, varying from only the single central point (i,j) for 
the point-Jacobi scheme to all five points for the ADI scheme. 

Convergence-rate estimates for each of these schemes are now presented. 
In brief, this estimation technique proceeds as follows (see Ames, ref. 9 for 
a more detailed discussion) . It can be shown that 

1I."'"‘II/IU'‘1I « » 

where is the error in the nth-level solution defined by 


and X is the average spectral radius of the iteration scheme (that is, the 
ratio of the average maximum eigenvalue at n + 1 to the average maximum 
eigenvalue at n) . To achieve a given reduction in the initial error, say 
p, we have 

p = Ile"ll/||e°ll a x" (3.46) 

Thus, -log^g p is. the number of orders of magnitude by which the error is 
reduced in n iterations. Rearranging, we have 

Nura = logjp X/log^u p (3.47) 

Equation (3.47) gives the number of iterations required to reduce the initial 
error, e° , by -log^^ p orders of magnitude. Ames (ref. 9) evaluates an 
expression similar to equation (3.47) for each of the relaxation schemes just 
presented. For a one-order magnitude drop in the error, that is, for 
-logji Q p = 1, an estimate of the number of iterations (Num) can be obtained 
and is given for each scheme in table 1 . The quantity A is the mesh spac- 
ing. Note that the line versions of each algorithm are about a factor of 2 
faster than the corresponding point counterparts (except SLOR which is 
faster than SOR) . More important, the SOR and SLOR schemes are much faster 
than the other point or line schemes. For example, when A = 0.01, SOR is 
200 times faster than point-Gauss-Seidel. Thus, the importance of overrelax- 
ation in relaxation schemes in emphasized. Furthermore, ADI is much faster 
than SLOR; for example, when A = 0.01, ADI is over 10 times faster than 
SLOR. These convergence rate estimates are compared graphically in figure 7. 
For more details regarding the assumptions of this analysis see Ames (ref. 9). 


(3.44) 

(3.45) 


4. ALGORITHMS FOR THE FULL-POTENTIAL EQUATION: EARLY IDEAS 


In this section, priliminary concepts associated with the numerical 
solution of the full-potential equation are discussed. Because the noncon- 
servative form of the full-potential equation is more revealing to linear 
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analysis and also because It was the first form of the full-potential equa- 
tion solved, we will restrict our attention solely to this equation in this 
section. The nonconservatlve full-potential equation was presented earlier 
(see sec. 2.3) and Is restated here for convenience: 

{a} - - 2uv(|i^y + (a^ - v^)(J>yy - 0 (4.1) 

For purposes of simplifying the analysis we can treat the coefficients 
a* - u^ , -2uv, and a^ - v^, as local constants. Thus, the equation Is 
effectively linearized. 


4.1 Subsonic Difference Scheme 

A spatial finite-difference scheme for the nonconservative full- 
potential equation (4.1), which is valid for subsonic flow only, is given by 


This scheme is centrally differenced and second-order accurate. If the 
point-Jacobi iteration scheme is now applied to equation (4.2) we have 




2(a' 


u^) 2(ai_- 






V^)1 pH 




(4.3) 


Stability can be investigated for this scheme by using the von Neumann 
test. The amplification factor (G) for Ax = Ay is given by 

G = 1 + R(cos aAx - 1) + S(cos aAx sin bAy) 

+ T(cos bAy - 1) (4.4) 

where 



and q is the speed of the fluid (/u^ + v^ ) • As before with Laplace's 
equation, the amplification factor is purely real. For stability, the mag- 
nitude of the amplification factor must be less than or equal to 1. When 
the flow is locally supersonic (q > a) , this centrally differenced scheme is 
unstable. An example demonstrating this situation is presented as follows: 
Assume the flow is aligned with the x-axis; therefore, u = q and v « 0. 
Equation (4.4) reduces to 

2 

G = 1 + R(cos aAx - 1) + (cos bAy - 1) (4.6) 

2a^ - q^ 

Since cos aAx - 1 < 0 and R < 0 for q > a, it is obvious that G > 1. 
That is, this centrally differenced scheme is unstable for supersonic flow. 

Proving stability for subsonic flow is a little more difficult. The 
amplification factor must satisfy -1 i G < 1. Working with only the first 
inequality, we have from equation (4.4) the following condition: 

R(1 - cos aAx) - S sin aAx sin bAy + T(1 - cos bAy) <2 (4.7) 
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The condition for subsonic flow requires the discriminate of equation (2,3) 
to be less than or equal to zero. This can be restated as 


✓St > |s 


(4.8) 


This expression is used to replace the S coefficient in the cross-product 
term and makes the original condition for stability more restrictive. The 
sign of this cross-product terra can be either positive or negative. Clearly, 
if this term is negative then the inequality is always satisfied. Therefore, 
assuming the sign to be positive yields 

R(1 - cos aAx) + v^r|sin aAx| [sla bAyj + T(1 - cos bAy) <2 (4.9) 

Next , add and subtract the terms necess» 9 iry to combine the cross-product term 
in a perfect square. This yields 

R(1 - cos aAx) + T(1 - cos bAy) + y sin^ aAx + ~ sin^ bAy 

-j2 

|sin aAx| - |sin bAy|j < 2 (4.10) 

Because the perfect-square term is always negative, it can be neglected, 
making the original expression more restrictive. After further simplifica- 
tion the inequality becomes 



- -2 (1 + cos aAx)‘ 


- y (1 + cos bAy)^ < 


(4.11) 


The inequality is now seen to always be true. Therefore, the first half of 
the stability analysis is complete. Looking at the second half yields 

R(cos aAx - 1) + S sin aAx sin bAy + T(cos bAy - 1) < 0 (4.12) 

Using a similar substitution on this expression to reduce the cross-product 
term yields, after simplification, 

- y (1 - cos aAx)^ - y (1 - cos bAy)^ ^ 0 (4.13) 

This inequality is also always true. Thus, the simple central-difference 
scheme given by equation (4.2) is stable, providing q < a, but is unstable 
when q > a. In the example just considered the iteration scheme was point- 
Jacobi; however, these conclusions can (theoretically) be reached with any 
iteration scheme. 


4.2 Simple Supersonic Difference Scheme 

A scheme suitable for supersonic regions of flow was first introduced 
by Murman and Cole (ref. 24), for the TSD equation. This approach was later 
advanced by many researchers including Steger and Lomax (ref. 25), Garabedian 
and Korn (ref. 26), Bailey and Steger (ref. 27), and Eallhaus and Bailey 
(ref. 28) for a variety of formulations in both two and three dimensions. 

The basic idea is as follows: First, determine the local flow type at each 

grid point (either elliptic or hyperbolic), by central differencing the 
velocity potential. Then, at subsonic points, use the standard central- 
difference approximation (second-order accurate), and at supersonic points 
use difference formulas retarded in the upwind direction. This yields an 


overall differencing scheme In which the computational domain of dependence 
correctly reflects the physical domain of dependence. 

When the flow Is essentially aligned with the positive x-coordlnate 
direction the following simple upwind difference scheme Is ultable for 
supersonic regions of flow (see Jameson, ref. 29): 



[(a^-u^)E“^6 - 2uvt 6 + (a^ 

X XX X y 




(4.14) 


Upwind evaluation of the term, as in equation (4.14), models the origi- 

nal second-order central-difference approximation plus an upwind-differenced 
artificial viscosity term. This can be shovm as 


E"^6 <t> - (6 - <S + E~^6 )d> 

X XX ' XX XX X xx“^^ 


■( 


1 - E 


-1 


6 )(fi 
X - X.J 


6 - Ax — 

XX Ax XX 


= (<S - Ax^ 6 )ij) 

XX X xx'’^ 


= <t>xx - AX.I.XXX + O(Ax^) (4.15) 

where the -Axiti^xx is the upwind-evaluated artificial viscosity. [Equa- 
tion (4.15) can also be derived by a standard Taylor series approach.] The 
spatial difference scheme given by equation (4.14) is only first-order accu- 
rate. Second-order accuracy In subsonic regions and first-order accuracy in 
supersonic regions are typical characteristics of most successful transonic 
difference schemes based on potential formulations. 


An appropriate iteration algorithm for solving the spatial difference 
scheme expressed by equation (4.14) is given simply by 


[(a^ - u^)E 6 

X XX 


2uvt 6. + (a^^ - v2)6_Jc" . + Lf/ . = 0 

» J 


X y 


yy 


*^i.j 


(4.16) 


where the N-operator used here is the same as the L-operator defined in 
equation (4.14). The velocity components used in both .the N and L oper- 
ators of equation (4.16) are defined by 


n 

u. . 


= 6 V . 

x^i.j 



6 d)' 

\ T ~ 


y^i.j 


(4.17) 


and a is obtained from these velocity components and equation (2.18). The 
iteration scheme of equation (4.16) requires the inversion of a set of tri- 
diagonal matrix equations along the y-direction and is most judiciously 
implemented by sweeping downstream along the x-direction. 

Numerical stability of the simple supersonic algorithm just presented 
can be investigated using the von Neumann stability analysis. Because the 
N and L operators are identical, the amplification factor is zero, indicat- 
ing a direct procedure which converges in one iteration. Even though the 
solution procedure in che supersonic region is very fast, it is not direct. 
There are two reasons for this. First, the equation being solved is actually 
nonlinear and must be iterated to convergence. This factor is not measured 
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by the linear stability analysis. Second, the sonic line boundary supplied 
to this scheme as initial data must be iterated until the proper position is 
obtained, a process which may take many iterations. Note that if an over- 
relaxation factor (m > 1) had been used for this supersonic marching scheme 
the convergence rate would have been slower, that is, the amplification fac- 
tor would not have been zero. Setting all relaxation factors to 1 for super- 
sonic regions is a standard procedure for most line relaxation schemes. 

The Iteration scheme given by equation (4.16) is essentially a hyper- 
bolic marching Scheme. Therefore, it is appropriate to check the marching 
stability of the scheme. This is accomplished by the usual von Neumann 
stability test in which 



ax iby 
e e 


(4.18) 


is substituted into equation (4.16) to yield, after simplification, a quad- 
ratic expression in Solving this expression for yields 


-aAx , 
e =1 


1 ^ sin bAy ± 




sin^ bAy - 2 -r (cos bAy - 1) (4.19) 


The A, B, and C coefficients are given by 



B 


uv 

Ax Ay ’ 



For stability 


G = 


aAx 


I i 1 


or equivalently 


-aAx 


e 


> 1 


(4.20) 


(4.21) 


(4.22) 


Clearly, the condition for stability will be satisfied if the quantity under 
the radical in equation (4.19) is always negative. The first term is always 
negative and the second term is negative providing C/A < 0. The quantity 
C/A is negative as long as the x-marching direction remains supersonic, 
that is, as long as u > a. (Keep in mind that this scheme has been formu- 
lated for flows essentially aligned with the x-axia where v is small with 
respect to u.) Thus, the supersonic marching scheme given by equa- 
tion (4.16) is stable, providing the x-direction remains hyperbolic. When 
the X marching direction becomes subsonic (u < a) even though the flow is 
still supersonic (q > a) , a marching instability is predicted. 

A schematic example of this situation, using characteristics, is shown 
in figure. 8. Notice that the computational domain of dependence does not 
include the entire physical domain of dependence, a situation which roust lead 
to trouble. Near sonic lines, when the flow is only slightly misaligned with 
the x-axis, instability is predicted. Theoretically, this makes the simple 
supersonic difference scheme given by equation (4.14) impractical. 

The marching instability can also be attributed to a sign change of the 
artificial viscosity term, Ax(u^ - a^)<|ixxx» which should be positive for 
stable operation. This term, as mentioned previously, iJ a consequence of 
the upwind bias on the differencing scheme in supersonic regions and is the 
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mechanism for allowing physically correct compr^*sslon shocks while disallow- 
ing the physically Incorrect expansion shocks < 

In practice, the simple supersonic difference scheme Is suitable for 
small regions of supersonic flow In which the flow Is nearly aligned In the 
proper direction. Small regions of flow near the sonic line with negative 
artificial viscosity are not enough to prevent stable operation (ref. 29). 
However, when the free stream Is supersonic or when a swept-wlng calculation 
la Involved, a new technique for handling supersonic regions is necessary. 


4.3 Rotated Supersonic Difference Scheme 

The concept of rotated differencing was first introduced by Jameson 
(ref. 29) in 1974 and used for solving the nonconservative full-potential 
equation for transonic flows about wing geometries. This concept utilizes 
the Ideas already discussed, namely, second-order-accurate central differenc- 
ing for subsonic flow and first-order-accurate upwind differencing for super- 
sonic flow. However, to remove the directional difficulties associated with 
the previous simple supersonic differencing scheme, the following coordinate- 
invariant difference scheme is used for supersonic regions. For two- 
dimensional cases the nonconservative full-potential equation can be trans- 
formed into a local stream and stream-normal coordinate system (s - n) , 
which yields 

“ 0 (4.23) 

the X and y Cartesian coor- 


(4.24) 


(4.25) 

ijigg and are expressed in terms of Cartesian-coordinate derivatives 

by 

(4.26) 

♦nn ■ ^ <*•”) 

[Note that the characteristic directions associated with eq. (4.23) are always 
symmetric about the stream direction; see eq. (2.3).] 

The algorithm proposed by Jameson (ref. 29), which is based on equa- 
tion (4.23), is now presented. At hyperbolic points (q > a), central- 
difference formulas are used for all contributions to [eq. (4.27)]. 

These formulas are given by 


(a^ - q^)^3s 


where the s and n coordinates are relattitd to 
dlnates by 
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X “ - s - - n 
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Equations (4.28) and (4.29) have been written assuming u > 0 and v > 0. 

If these signs are different, equations (4.28) and (4.29) are replaced by 
similar formulas which retard the difference scheme in the proper upwind 
direction. The n-superscripts have been determined so as to provide favor- 
able temporal damping to the iteration scheme in supersonic regions. More 
discussion on this topic will be presented shortly. 

( 

The iteration scheme obtained by substituting equations (^f. 28) and (4.29) 
into equation (4.23) is expressed in standard correction form by 
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(4.30) 


- * ....W Ygg 

differencing, and the terms inside the second set of square brackets arise 
from the differencing. The spatial differencing molecules associated 

with this scheme are shown in figure 9. All four variations corresponding to 
the velocity vector originating in the four quadrants are shown. The contri- 
bution due to the term is shown by + symbols , and the contribution 

owing to the i|)gg term is shown by o symbols. Note that the <jinn differ- 
encing molecule never changes but that the 4isa differencing molecule is 
automatically upwind biased into the quadrant from which the velocity arises. 
The transition between these differencing variations is smooth because as one 
difference operator is switched to another the lead coefficient of that term 
automatically goes through zero. 

The rotated difference scheme just presented, completely removes the 
marcl^ktig instability problems associated with the simple upwind scheme pre- 
sented earlier. However, several new limitations or disadvantages are intro- 
duced: (1) the rotated differencing scheme no longer mimics a direct march- 

ing scheme and, therefore, slows convergence somewhat; (2) the increased size 
of the computational molecule, as well as the fact that the scheme is first- 
order accurate in both directions, increases shock smearing; and (3) this 
scheme must be swept in the flow direction, a minor limitation associated 
with general curvilinear meshes. Despite these limitations, this scheme 
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pioneered the way for three-dimensional, transonic flow-field calculations 
using the full-potential equation and is still the basis for a widely used 
computer code in the aircraft Industry — F1022. 


4.4 Temporal Damping 

The equivalent time-dependent equation associated with the iteration 
scheme of equation (4.30) can be written as 


2a<(.^^ + + (M^ 


l)(j> - (j, - 0 

' ^88 ^nn 


(4.31) 


where the a and @ coefficients are defined by 


a - (M^ - 1) M ^ V 
Ax q Ay q 

p _ 1 At V 

2 Ax q 


(4.32) 

(4.33) 


The quantity, (jigt, built into the iteration scheme of equation (4.30), pro- 
vides temporal damping. This is very Important in supersonic regions and is 
required to maintain stability. The term must be differenced upwind 

and with the proper sign so as to add to the magnitude of the matrix diagonal. 
This ensures diagonal dominance for the matrix inversion. 

Introducing an independent -variable transformation 

t » T + aX/(M* - 1) - 0Y 


s = X (4.34) 

n = Y 


for equation (4.31) yields 




2 



+ (M' - 


- 


0 


(4.35) 


For supersonic flows, it can be seen that the coordinate T is space-like 
(elliptic) and that either X (which is the streamwise coordinate) or Y 
(which is the stream-normal coordinate) is time-like (hyperbolic). In the 
physical steady-state problem, X is time-like; therefore, to maintain this 
situation computationally, the following condition is required 

> B^(m2 - 1) (4.36) 


Rewriting equation (4.36) using the definitions of a and B [eqs. (4.32) 
and (4.33)1 yields 


/At u 

^ At V 

>1 

/At 

\Ax q 

Ay q/ 

4 

VAx 


(4.37) 


Generally speaking, the iteration scheme of equation (4.30) has been designed 
to satisfy the condition of equation (4.37). However, problems arise near 
sonic lines (as - 1 ->■ 0). Also, in many iteration schemes the <^gj. tem- 

poral damping term is not automatically added. For both of these cases (Jigj. 
can be added explicitly. A suitable term for the present case is given by 
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(4.38) 


where the <)ixt '^yt terms are differenced as follows (u > 0, v > 0): 



^ 1 


.n . .n+i 

- h.i + ♦i-i.j 

A** 

~ *i-i,j 

xt 

Atilx 


„ 1 

^*i.j 

• n , ,n+i 

" '*’l.j ^ ‘*’i.j-i 

a" 

'yt 

AtAy 


(4.39) 

(4.40) 


The parameter e in expression (4.38) Is a user-specified constant, set 
large enough for stability but not so large as to excessively slow 
convergence . 


5. TRANSFORMATION AND GRID-GENERATION TECHNIQUES 


5.1 Governing Equation Transformation 

So far, all examples presented have dealt with only Cartesian coordi- 
nates. Usually, before solution algorithms can be implemented, the governing 
equation must be transformed from Cartesian coordinates into some suitable 
computational domain. Even applications that use Cartesian coordinates in 
the computational domain typically require the use of stretching or shearing 
transformations or both. The primary reason for applying an Independent 
variable transformation to the governing equation is to transform any geomet- 
rical surfaces in the problem into constant coordinate lines in the computa- 
tional domain. Thus, boundary-comdition implementation and grid clustering 
at geometrical surfaces can be achieved without undue difficulty. 


A general, independent variable transformation typically used in con- 
junction with numerical grid-generation procedures, is given by 


5 “ ?(x,y) 
n = n(x,y) 


(5.1) 


where x and y represent the Cartesian-coordinate physical domain and C and 
T) the computational domain (see fig. 10) . The conservative full-potential 
equation written in Cartesian coordinates was presented in section 2 and is 
restated here for convenience as follows: 


'*■ ® (5.2a) 

p = I' - (<^x + ^5>] <5-2b) 

Equation (5.2) transformed by equation (5.1) yields the full-potential equa- 
tion written in the (?,ri) computational domain. 


P 


(f \ ^ (f )„ ■ 0 

[l . ^ ' 


(5.3a) 


(5.3b) 
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where 


and 


U 

V 





(5.4a) 


(5.4b) 

VC • 

(5.5a) 

VC • Vn - Cj^n^ + KyUy 

(5.5b) 

Vn • Vn - + n^ 

* y 

(5.5c) 


? n - ? n 

X y y X 


(5.6) 


The quantities U and V are contravariant velocity components along the 
5 and n directions, respectively; Aj, A2 , and A3 are metric quantities; and 
J is the Jacobian of the treq- '^tion. To evaluate the metric quantities 

of equations ( 5 . 5 ) and ( 5 . 6 ), '5 fallowing metric identities are required: 


J = l/(x^y^ - x^y^) 

?x = “ -Jy? 


(5.7) 


( 5 . 8 ) 


This transformation maintains the strong conservation-law form of the 
original equation and hence possesses characteristics suitable for a shock- 
capturing scheme. Additional Information about such transformations can be 
found in references 30 - 34 . 


The transformation metric quantities just presented — A^, A2 , A3, and 
J — all have interesting physical interpretations. First of all, the 
Jacobian, J, can be shown to approximate the inverse of the cell area (in 
three dimensions the Jacobian approximates the inverse of the cell volume) . 
Thus, a simple test on the sign of the numerically computed Jacobian will 
reveal if the mesh crosses over on itself (an occasional symptom of some 
numerically generated meshes which is devastating to the stability of the 
iteration scheme) . 

4 

In general, the (x,y) ■> (?,n) transformation produces a nonorthogonal 
mesh in the physical domain. A direct measure of the amount of mesh skewness 
or nonorthogonality is the A2 metric. Since A2 = V? • Vn, the special case 
of an orthogonal mesh is realized when all the A2 metrics are zero. The 
A^ and A3 metrics provide a measure of cell aspect ratio. The Aj^ metric 
is approximately the ratio of arc length along the n-direction to the arc 
length along the ^-direction. The A3 metric is approximately the Inverse 
of A^. Thus, intimate details regarding a finite-difference mesh can be 
obtained automatically, without even plotting the grid, by Just monitoring 
the A3, A2 , A3, and J quantities. 

The three-dimensional form of the full-potential equation written in 
general curvilinear coordinates (?,n,?) is given by 
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The contravarlant velocity components U, V, and W are defined by 


(5.9a) 


(5.9b) 


(U,V,W)^ - h'“^H(«^^,(|.^,(!,^)'^ 


the quantity H Is defined by 


, j.(,€.»il„v?) 
8(x,y,z) 




’^x \ \ 


y zy 


(5.10) 


(5.11) 


and J » det(H). In equation (5.11), x, y, and z are the usual Cartesian 
coordinates along the flow, span, and vertical directions, respectively. The 
Til and ^ computational coordinates represent wraparound, spanwlse, and 
radlal-llke directions, respectively. 


5.2 The Elliptic Grid-Generation Procedure 


Perhaps the most popular new technique for grid generation Is the numeri- 
cal approach. The main idea associated with this technique is to establish a 
set of curvilinear coordinates by requiring that they be solutions to appro- 
priately formulated PDEs. The properties of the PDEs are such that srol»th 
and regular finite-difference grids result. Coordinate lines of the 
family do not cross, and coordinate linen of opposite families int'^ ’ 
orthogonally or nearly orthogonally. Dirichlet boundary conditions 4 ?:e speci- 
fied such that the body or other desired boundaries are automatically mapped 
to constant coordinate lines in the physical domain. This guarantees a one- 
to-one mapping in which the mesh is well-ordered and body-fitted. Mechanisms 
for achieving general mesh control exist through boundary-condition specifi- 
cation and through the control of various arbitrary coefficients depending on 
exactly which formulation is used. The numerical mapping procedure is 
generally valid for both two- and three-dimensional flows. 


Many numerical grid-generation techniques have been developed for a wide 
variety of applications. However, only three schemes will be discussed here. 
For an extensive survey of this subject see Thompson, et al. (ref. 35). The 
first scheme has its roots in a number of publications, but has been recently 
developed and popularized by Thompson et al. (refs. 36-38). This scheme is 
perhaps the most widely used numerical grid-generation scheme. Poisson's 
equation is used to define the transformation and is given by 


?xx + Sy = 

\x \y = Q<5,n)) 


(5.12) 


The P and Q right-hand-side quantities are defined as a sum of exponential 
terms each with several free coefficients. These coefficients can be adjusted 
by the user to provide different types of mesh size and skewness control. 
Equations (5.12) are transformed to (and numerically solved in) the 
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computational domain, that Is, the roles of (C,n) and (x,y) are Interchanged 
such that equations (5.12) become 






ay 




^nn “ 


+ Qx^) I 

+ Qyn) I 


(5.13) 


where 


2 2 
a =» X + y 
n ■'n 




Xc-y« “ XV- 

? n n ? 


4* 


(5.14) 

(5.15) 


The finite-difference grid is generated by numerically solving equa- 
tions (5.13)-(5.15) . First, all derivatives are replaced by standard second- 
order-accurate finite differences. The spatial Increments, A? and An, can be 
arbitrarily chosen and are usually set to 1. Once boundary and Initial 
values of x and y are specified, the final interior values can be computed 
by standard relaxation procedures. 


An important aspect associated with numerical elliptic-solver techniques 
is that they have a high degree of controllability. The vast number of free 
coefficients contained in the P and Q terms is an indication of the large 
amount of control available to the user. There is an obvious difficulty 
associated with this flexibility; How can it, in a general and simple way, 
be made available to the user? This problem poses a difficulty in two 
dimensions and becomes seemingly insurmountable in three dimensions. 


Several researchers have experimented with different aspects of this 
problem, including Thomas and Middlecoff (ref. 39), Steger and Sorenson 
(ref. 40), and Sorenson and Steger (ref. 41). In the latter approach, a 
simplified form of the P and Q terms was adopted and is given by 


p - Pocor*'"-"') 


where corresponds to the hmin inner boundary, Pg and Qo are sets of 

constants which vary in the ? direction on the hniin boundary, and a and 
b are constants which control the rate of decay of the P and Q forcing 
terms into the mesh interior. Thus, application of the forcing terras is 
restricted only to where the control is primarily desired. Steger and 

Sorenson recognized that two types of control are desired: (1) control of 

the cell aspect ratio at boundaries and (2) control of cell skewness at 
boundaries. These two conditions can be expressed mathematically by 


s = 


ds 

dn 


(< 


+ y"") 


1/2 


As 


(5.17) 


and 

• Vn = x^x^ + y^y^ = |V5| |Vn|cos 0 (5.18) 


where s is the arc length along the q lines which Intersect the body and 
9 is the angle of intersection., The condition for orthogonality Is 6 * ir/2. 
Other angle specifications are possible, thus allowing precise control on the 
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cell skewness at the grid boundary. By utilizing equations (5. 16)-(5. 18) , 
as well as the Dirichlet boundary conditions in x and y, expressions for 
Pq and Qq can be derived. Specification of As(?) and 0(§) is all that is 
required to completely determine a new grid with the desired cell aspect 
ratio and cell skewness at the boundary. 

An example of this technique used to generate a grid about a highly cam- 
bered 12-to-l ellipse is shown in figures 11 and 12. Figure 11 shows global 
and trailing edge detail of the Laplacian mesh (that is, a mesh with no con- 
trol) and figure 12 shows the same views with control (9 = tt/2 and 
As ■ 0.005). Note the poor grid quality for the case with no control, espe- 
cially in the concave portion of the ellipse. The case with control, how- 
ever, produces a nearly orthogonal grid at the body with a nearly uniform 
As distribution closely approximating the desired value of 0.005. 


5.3 The Hyperbolic Grid-Generation Procedure 

The next numerical grid-mapping procedure to be discussed is based on a 
hyperbolic set of governing equations. This type of grid-generation proce- 
dure is less well developed than the elliptic procedures but does have sev- 
eral desirable characteristics. Examples of this type of grid-generation 
procedure include Starius (ref. 42) and Steger and Chaussee (ref. 43). In 
the latter work, the PDEs used to numerically define the finite-difference 
mesh are given by 


vc*vn = 5n + ?n = <t> 
XX y y 


? n - ? n = J 
X y y X 


(5.19) 

(5.20) 


Orthogonality is achieved by simply setting (p to zero, and J represents 
the Jacobian of the transformation (that is, effectively the grid-cell area). 
In this formulation, both p and J are user-specified functions. This pro- 
vides a great deal of controllability which is perhaps more direct and easier 
to implement than in most elliptic PDE formulations. Direct specification of 
J produces a well-behaved mesh that does not cross over on Itself except in 
the most severe cases. Because this system of equations is hyperbolic in q, 
a solution can be obtained simply by marching away from initial data speci- 
fied on the inner boundary. Because iteration is not required, computation 
time for this technique is very small. 


A disadvantage of this technique lies in the lack of direct control over 
the position and distribution of grid points on the outer boundary (n = n nia-w ) 
However, for external aerodynamic applications this limitation is not severe. 
Other more fundamental problems lie in the treatment of surface singularities 
or extension to three-dimensions. Because of the hyperbolic nature of the 
grid-transformation equations, any singularities imposed in the n = Bmin 
initial data — for example, a slope discontinuity in the geometry — will 
propagate in the n-direction and perhaps cause problems in the grid interior 
Formulations of a suitable algorithm in three dimensions is more difficult 
because the orthogonality condition [eq. (5.19)] expands to three equations. 
Thus , a system of four independent equations is created where only three can 
be used. 


Much work needs to be completed before numerical grid-generation schemes 
based on hyperbolic PDEs can be used routinely for three-dimensional calcula- 
tions. However, with the advantages of speed and controllability, this grid- 
generation procedure is an attractive alternative to the elliptic-solver 
grid-generation technique. 
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I 5.4 The Parabolic Grid-Generation Procedure 

I The final numerical grid-generation procedure presented in these notes 


^ la now discussed. This procedure, developed by Nakamura (refs. 44 and 45), 

, is based on a parabolic set of equations given by 

i 

^ A tS'ts-. ■ ^ n « 


1 

X » Ax,-c - 2Bx- + S 
n Cn X 

(5.21a) 


y » Ay__ - 2By^ + S 
•'n ■'e? Cn y 

(5.21b) 

! . ■ 

1 

where A and B are positive constants, ^ and n are the usual computational 
coordinates along and away from the airfoil surface, x and y are the physi- 
cal domain Cartesian coordinates, and and Sy are source terms given by 

i 

- 3) 

(5.22a) 


^ll.J ■ ('’l.HJ - - 3) 

(5.22b) 


In equations (5.22), the i and j subscripts represent position in the 
finite-difference mesh (C = iAC and n = jAn) and j = NJ is the outer 
boundary. Thus, the effect of the outer boundary (x^^jjj and yi^Nj) 
included in the Sj^ and Sy source terms. 

Equations (5.21) are parabolic in n and can be marched away from ini- 
tial data (xi I and y^ i). These initial coordinates represent the airfoil 
surface, specified by the user, with any smooth distribution. A standard 
i difference scheme involving centered differences for all 5 derivatives and 

i backward differences for all n derivatives is used to discretize equa- 

i tions (5.21). This produces a set of uncoupled tridiagonal matrix equations, 

I which, when inverted, yield values of x and y at j * 2. This process is 

repeated marching away from the airfoil surface, j =1, to the outer boundary, 
* j = NJ. As j approaches NJ the x and y values automatically approach 

j the specified outer boundary values of xi^jjj and yi,NJ* 

I An example of this procedure applied to a highly cambered thin ellipse 

I is shown in figure 13 (taken from ref. 44). A view of the entire geometry is 

i shown in figure 13(a) and a close-up of the trailing edge is shown in fig- 

' ure 13(b). Note the regularity of cell size around the inner boundary and 

; lack of cell skewness over the entire mesh. This grid compares quite favor- 

j , ably with the grid presented in figure 12 which was generated with explicit 

I controls on cell size and skewness by the elliptic-solver approach. However, 

I the parabolic scheme requires only a fraction of the computational work. 

I 

I The parabolic scheme is also easily extended to three dimensions (see 

I ref. 45). One additional equation is required to define the third coordinate 

\ and each equation has additional terms; otherwise, the extension is quite 

straightforward. Because the resulting difference equations are no longer 
narrow-banded an ADI factorization is used. Thus, in proceeding from one 
coordinate surface to the next, several iterations are required. However, 
the basic parabolic grid-generation scheme is still direct, requiring only 
one sweep from the inner to the outer boundary and, thus, requires only a 
fraction of the computer time used by the elliptic-solver scheme. 

i An example of a three-dimensional grid generated by the parabolic scheme 

about a wing/fuselage geometry is shown in figure 14 (taken from ref. 45). 

The wing planform plane, the fuselage surface, and a typical fuselage cross- 
sectional grid are all highlighted in this figure. 
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5.5 Geometric and Analytic Mapping Procedures 

Algebraic and analytical mapping procedures can be simple ~ for example, 
Involving only a stretching or shearing of the coordinate system. On the 
other hand, they can be more complicated — for example, conformal mappings. 

The stretching and shearing transformations can be useful for simple geome- 
tries but generally are not sufficient by themselves for more complex airfoil 
and wing calculations. Conformal mappings do have the generality required 
for providing good quality, economical grids for reasonably complex, two- 
dimensional geometries. Many researchers have used conformal mappings to 
generate arbitrary, orthogonal (or nearly orthogonal meshes when a sheared 
conformal mapping procedure is used) for a host of two-dimensional applica- 
tions. A few examples include those of Sells (ref. 46) for an airfoil; 
Kacprzynski (ref. 47) for an airfoil between wind-tunnel walls; Ives (ref. 48) 
for a multielement airfoil; Caughey and Jameson (ref. 49) and Chen and Caughey 
(ref. 50) for axisymmetric Inlets with and without a centerbody; and Ives and 
Llutermoza (ref. 51) for axial-flow turbomachinery cascade applications. 

The theory behind conformal mapping techniques is governed by analytic 
functions of a single complex variable. This theory is well developed, but 
fundamentally limited to applications in two space dimensions. Nevertheless, 
some researchers Iv.ve found ways to use conformal mapping techniques to 
assist in generating grids for three-dimensional problems. A few examples 
are given by Jameson (ref. 19) for wings; Jameson and Caughey (ref. 52) and 
Caughey and Jameson (refs. 53 and 54) for wind/body combinations/ and Ives 
and Manor (ref. 55) for three-dimensional inlet and inlet-centerbody config- 
urations. The basic approach utilized for three-dimensional grids is to 
generate a series of two-dimensional grids using standard conformal mapping 
procedures. Then these grids are "stacked" together in the third dimension 
to form the final three-dimensional grid. This approach has worked well for 
geometries with smooth variation in the third dimension but lacks generality. 
Conformal mapping procedures will continue to be used successfully in a wide 
variety of applications; however, the anticipated trend for grid-generation 
procedures will be away from such techniques and toward the more general 
(although presently less well-developed) numerical or geometric mapping 
procedures . 

The geometric grid-generation procedure has recently received much 
attention and promises to develop into a practical approach. These proce- 
dures are efficient, requiring very little computer time, and have general 
capabilities regarding coordinate line control. In addition, this type of 
procedure seems to be general enough for easy extension to three dimensions. 
Eiseman has presented geometric grid-generation techniques for several two- 
and three-dimensional configurations (refs. 56-58). The mathematical aspects 
of this new geometrical grid-generation procedure are developed by Eiseman in 
references 59 and 60. An example application using this grid-generation pro- 
cedure for the numerical computation of transonic airfoil flows is given by 
Pulliam et al. (ref. 61). 

Other geometric grid-generation procedures have been presented by 
McNally (ref. 62), Graves (ref. 63), and Eriksson (refs. 64 and 65). The 
latter procedure, which can be viewed as a generalized spline interpolation 
technique, generates a finite-difference grid in two or three dimensions by 
interpolating geometric data from the domain boundaries. The geometric data 
required include boundary coordinates and derivatives. In some cases outer 
domain boundaries can automatically be determined by the geometric transfor- 
mation procedure and, therefore, do not necessarily have to be specified. 

This procedure is computationally efficient and has good coordinate line con- 
trol properties. Transonic wing computations using the Euler equations have 
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been obtained using this grid mapping procedure for a variety of configura- 
tions (ref. 66). 

The geometric grid-generation procedure has only recently been intro- 
duced and, therefore, its ultimate role is difficult to predict. However, 
the general properties of this class of grid-mapping procedures (e.g., flexi- 
bility, computational efficiency, and controllability) suggest a large 
increase in its use for the generation of grids about complicated 
configurations . 


5.6 Solutions-Adaptive Grid-Generation Schemes 

A solution-adaptive grid (SAG) technique is defined to be a grid- 
generation procedure in which the flow-field solution influences the genera- 
tion of the grid. Usually this solution influence is designed to produce 
mesh-point clustering about flow-field gradients and, therefore, produce 
reduced levels of truncation error. Several researchers have experimented 
with different types of SAG algorithms, including Dwyer et al. (ref. 67), who 
developed a procedure for various time-accurate heat-transfer problems; 
Glowinski (ref. 68), who experimented with optimal grids for incompressible, 
inviscid flow, using a finite-element technique; and Pierson and Kutler 
(ref. 69), who determined optimal grid-point distributions for several model 
problems. Additional work has been presented by Rai and Anderson (refs. 70 
to 72), in which the positions of the grid points are determined iteratively 
along with the solution-dependent variables. The method used to adjust the 
grid-point positions , forces the solution gradients to be averaged or spread 
evenly over the finite-difference mesh in the computational domain. This has 
the effect of reducing the average truncation error. 

Other examples of SAG approaches applied to transonic airfoil solutions 
using the full-potential equation are presented by Ushimaru (ref. 73), Holst 
and Brown (ref. 74), and Nakamura and Holst (ref. 75). In the latter 
approach, the surface grid distribution is redistributed by using a second- 
order ordinary differential equation (ODE) in which the dependent variable is 
a grid-density function. The forcing function for this ODE is a function of 
the airfoil surface solution gradient which has to be supplied from an inde- 
pendent solution computed on a standard mesh. In regions of large gradient 
(for example, at shock waves) the surface grid distribution is automatically 
clustered, and in regions of small gradient the surface grid distribution is 
coarsened. The interior grid is generated via interpolation using character- 
istics of the airfoil surface distribution to appropriately cluster the mesh 
interior. 

An example SAG grid taken from reference 75 is shown in figure 15(a). 
This grid contains 99 x 25 points and was generated about an NACA 0012 air- 
foil at a free-stream Mach number of 0.75 and an angle of attack of 2°. For 
convenience, only the inner 14 lines are shown. The surface-pressure coeffi- 
cient distribution computed with this grid is compared with a surface solu- 
tion generated on a much finer standard grid (245 x 56) in figure 15(b). 

Even though the fine standard grid contains 5.5 times more grid points, the 
shock is slightly steeper for the SAG solution. For this case despite the 
fact that the SAG procedure utilized two complete solutions, one on the stan- 
dard grid and the second on the SAG grid, the SAG procedure was 2 to 3 times 
faster than the fine-grid calculation. 
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6. SPATIAL DIFFERENCING SCHEMES 


6.1 The Finite-Volume Scheme 


The finite-volume spatial difference scheme of Jameson and Caughey 
(refs. 52-54; 76-77) was first used to solve the conservative full-potential 
equation in 1977. Since then many applications of this scheme have been 
made in both two and three dimensions. This scheme, written in two dimen- 
sions (for convenience) is given by (see ref. 77) 



i-bl /2 , j +1 


/2 


tt (^) 

^ 4+1/2, j + 1 /2 


= 0 


( 6 . 1 ) 


where the averaging (y^ and y_) and differencing (^^ and ^i^) operators are 
defined in section 3.1. The derivatives of x, y, and with respect to 
C and n are required to compute the density (p), the contravariant velocity 
components (U and V), and the Jacobian (J). These computations are all per- 
formed at cell centers (i+l/2,j+l/2) by using 


^ 44 + 1 / 2 , 3 + 1/2 4+1, 3+1 

^ 44 + 1 / 2,3 + 1/2 “ 4+1, j + 1 


(6.2a) 

(6.2b) 


The resulting spatial differencing scheme is very compact and requires only a 
single density evaluation per grid point. 


However, this scheme has a tendency to produce oscillatory solutions in 
which the i + j odd points are decoupled from the i + j even points. 

This situation can be corrected by adding suitable recoupling terms. The 
resulting scheme becomes 


n cVj 4 + 1 / 2 , j+ 1/2 ^ 4 + 1 / 2 , j+ 1/2 




(6.3) 


where 


4+1/2, j+1/2 J 

A I = — (A. 

n I i+i /2 , j+1/2 J ' 3 


j+1/2 

V^/a^)^+i/2,j+j^/2 


(6.4) 


A value of one half is generally used for the constant e. Addition of these 
terms recouples the odd and even points and represents a suitable spatial 
differencing scheme for subsonic regions of flow. 


The finite-volume scheme is stabilized in supersonic regions by the 
explicit addition of artificial viscosity terms given by 



(6.5a) 


(6.5b) 

where the switching function a is defined by 

0 = max[0, 1 - (M^/M)^] 

(6.6) 
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The Mq parameter used In equation (6.6) is a critical Mach number, defined 
in such a way that the spatial differencing scheme uses the subsonic differ- 
encing for values of local Mach number below M^ and the supersonic differ- 
encing for values of the local Mach number above In other words, the 

transition from central to upwind differencing does not necessarily take 
place at the sonic line. Note that M^ must be less than or equal to 1 for 
stability. 


The final spatial differencing scheme that is valid for both subsonic 
and supersonic regions of flow with the odd-even error compensating terms 
Included is given by 
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where the P and Q terms are defined as follows: 

'P 


■ i+1/2, j 


^i,j+i/2 
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i+i,j 


r^i,j+i 


if U > 0 
if U < 0 
if V > 0 
if V < 0 


(6.8a) 


(6.8b) 


The spatial differencing scheme given by equations (6.5)-(6.8) is cen- 
trally differenced and second-order accurate in subsonic regions where 
M < M(, and is upwind-differenced and first-order accurate in supersonic 
regions (or in regions where M > M{,). The upwind influence is retained for 
general curvilinear meshes regardless of the orientation of the velocity 
vector. Therefore, this conservative spatial differencing scheme approxi- 
mates the rotated differencing scheme first developed by Jameson (ref. 29) 
for the nonconservative form of the full-potential equation. 


Extensions of this scheme to higher orders of accuracy have been investi- 
gated by several researchers, including the work reported in references 76 
and 78-81. For example, Caughey and Jameson (ref. 76) modified the P and Q 
terms such that (looking at only the P term) 


i+1/2 , j 


Pi_j - (1 


if 


- ^V^^i+2.j 


U i 0 
U < 0 


(6.9) 


where c is a constant of order unity. When the local solution is smooth, 
the overall scheme is formally second-order accurate. In regions of high 
solution gradient (for example, at shock waves) these added terms force the 
overall accuracy to revert back to first order. This type of hybrid scheme 
has been found to be useful for maintaining stability in strong shock 
calculations. 
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6.2 The Artificial Density Scheme 

The artificial density spatial differencing scheme for the full- 
potential equation has been independently presented in several different 
forms (refs. 82-84). These forms, although not identical, have certain simi- 
larities which can be attributed to the earlier work of Jameson (ref. 85). 
Jameson's work is characterized by a scheme with an explicitly added artifi- 
cial viscosity term. This term biases the spatial difference scheme in the 
upwind direction for supersonic regions of flow but does not affect the 
centrally differenced scheme in subsonic regions. The three schemes of ref- 
erences 82-84 use this approach with one basic simplification: the upwind 

bias is accomplished by an upwind evaluation of the density coefficient. All 
three procedures compute this upwind or artificial density quantity in 
different ways. 


In the procedure of Holst and Ballhaus (ref. 83; see also refs. 34, 86) 
the finite-difference approximation for the full-potential equation written 
in general curvilinear coordinates [see eq. (5.3a)] is given by 



i+x /2 


.j 



,j + l/2 


0 


( 6 . 10 ) 


where the operators ^^( ) and 1i^r^( ) are first-order-accurate, backward- 
difference operators in the C and n directions, respectively (see sec. 3.1), 
and the density coefficients p and p are defined by 




(6.11a) 


Pi,j + l/2 “ “ '^^Pk,j+l/2 '^i,j + l/2Pi,j+S+l/2 


(6.11b) 


The r and s subscripts used in equations (6.11) control the upwind direc- 
tion of the density coefficients and are defined by 


r = ±1 when $ 0 

s = ±1 when ^ 0 


( 6 . 12 ) 


The switching or transition function V depends on the local Mach number 
Mi 4 and the flow direction and is defined by (e.g., looking at only the 
x-direction) 


'^i+ 1/2 ,j 


max[(Mj^^ - 1)C,0] for > 0 

Vi/2,j < 0 


(6.13) 


The quantity C is a user-specified constant usually set to a value between 
1 and 2. 


The density calculation is performed in a straightforward manner by using 
a discretized version of equation (5.3b). Values of the density are computed 
and stored at half points (i.e., at i + 1/2, j). Values of and 
required for computing the density at i + 1/2, j are given by 




(6.14) 
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The contravariant velocity components (U and V) used in equations (6.10) 

to ( 6 . 12 ) are computed with standard, second-order-accurate, finite-difference 

formulas. These quantities are computed and used at i + l/2,j or 

i»j + 1/2. An example differencing formula for U computed at 1 + 1/2, j 

is given by 

^i+l/ 2 ,j “ ■'^ili+i/2,j^5**’l,j ^2li+i/2,j^5^n'*’i.j 

The metric quantities A 2 , and A 3 are computed with fourth-order- 

accurate, finite-difference formulas and are stored at integer points in the 
finite-difference mesh. Values required at i + 1/2, j, i,j + 1/2, etc., are 
obtained by using simple second-order averages. This metric calculation pro- 
cedure has produced good results on smooth meshes but suffers on meshes that 
are not smooth. A superior metric differencing scheme which produces good 
results on even nonsmooth meshes will be discussed later in this section. 

With the spatial differencing scheme just outlined, an upwind Influence 
in supersonic regions is achieved without the explicit addition of an arti- 
ficial viscosity term. Instead, the stabilizing upwind influence is produced 
by the upwind evaluation of the density in an otherwise centrally differenced 
scheme. This approach is significant because it simplifies the technique for 
Including an upwind influence into the residual operator. As in the finite- 
volume scheme presented in section 6 . 1 , the present artificial density scheme 
closely approximates the effects of a rotated differencing scheme. This 
aspect greatly contributes to the stability and reliability of the present 
algorithm for many difficult test cases. 

Another variant of the artificial density spatial differencing scheme 
has been presented by Hafez et al. (ref. 84). In this scheme, which is desig- 
nated as an artificial compressibility scheme, the density coefficients in 
both coordinate directions are defined by 

"i.j ■ »i.j - '■i.j >i.j 

where 

> J » J 

The double-arrow notation indicates a first-order, upwind difference, s is 
the streamwise coordinate direction, and v is a switching function defined 
similarly to equation ( 6 . 6 ) or (6.13). 

One difficulty associated with the artificial density spatial discretiza- 
tion philosophy is that if either the switching function v or the density 
p (or both) are not properly computed, the shock capture process will produce 
large pre-shock oscillations and poor algorithm reliability. Two guidelines 
offered by South and Jameson (1979, private communication) that help elimi- 
nate this unacceptable behavior are ( 1 ) the quantity v should be evaluated 
at i,j not at i + l/ 2 ,j as in some formulations; and ( 2 ) the density 
values used in equation ( 6 . 11 ) or (6.16) should be computed at i+l/ 2 ,j+l/ 2 . 
The last guideline produces a density computation with a minimum-width dif- 
ferencing module in the streamwise direction. With these two guidelines the 
existence of pre-shock oscillations is greatly reduced. 

Examples showing the magnitude of these effects are shorn in fig- 
ures 16-18. Figure 16 shows two pressure coefficient distributions from 
reference 86 for an NACA 0012 airfoil at a free-stream Mach number of 0.75 
and an angle of attack of 2°. The two solutions correspond to different 


(6.16) 


(6.17) 
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values of C [see eq, (6.13)]. For both of th^se calculations, V Is incor- 
rectly evaluated at 1 + 1/2, j Instead of i,j, and p is computed 
directly at grid points instead of at the preferred location (i+ 1/2, j +1/2) . 
As a result a pre-shock oscillation exists even for relatively large values 
of C. Figure 17 shows the same comparison with v correctly evaluated at 
i,j. The pre-shock oscillation has been eliminated, even for relatively 
small values of C, but in some cases could still exist for values of v 
computed from equation (6.6) (when = 1.0). Figure 18 shows a pressure 
coefficient comparison from South et al. , (ref. 87) for a nonlifting NACA 
0012 airfoil calculation at a free-stream Mach number of 0.85. The two 
curves correspond to nodal-point (i,j) and mid-cell (f 4 1/2, j + 1/2) density 
calculations. The v parameter is correctly evaluated at i,j in both 
cases and defined by equation (6.6) with = 1.0. The mid-cell density 
calculation clearly gives the superior result with no oscillations. 

Many researchers have used one of the artificial-density, spatial- 
differencing schemes mentioned above because of the simple, reliable way in 
which the supersonic region is stabilized. A few of these applications 
include those of Farrell and Adamczyk (ref. 88), Akay and Ecer (ref. 89), and 
Deconinck and Hirsch (ref. 90), for cascade calculations; Shankar (ref. 91) 
for supersonic marching problems; Eberle (refs. 92 and 93) for a variety of 
different applications; and Steger and Caradonna (ref. 94) and Goorjian 
(ref. 95) for unsteady calculations. Results comparing a number of artifi- 
cial density scheme variations applied in a finite-element context are pre- 
sented in Habashi and Hafez (ref. 96). 


6.3 Spatial Differencing Schemes Based on Flux-Vector Splitting 

Spatial discretization schemes based on the flux-vector splitting models 
of Godunov (ref. 97) and Engquist and Osher (ref, 98) have recently gained in 
popularity. Their value in capturing shock waves sharply and in providing 
good stability properties for a variety of different iteration schemes has 
been reported by Goorjian and Van Buskirk (ref. 99), Goorjian et al. 

(ref. 100), Boerstoel (tef. 101), and Slooff (ref. 102). A comparison of 
these schemes is presented in van Leer (ref. 103) as they apply to the one- 
dimensional Burger's equation, 




F = 0 

X 


(6.18) 


where u is the flow velocity, and F is the flux, u^/2 (or for the full- 
potential equation F = pu) . These spatial discretization schemes applied to 
the steady part of equation (6.18) can be stated as follows (see ref. 102): 


where 


^J+l /2 = - "^^^4-i/2'^i+i/2> (Godunov) 

^1+1/2 = - 4-1/2 " 4+1/2 (Engquist-Osher) 


(6.19) 


( 6 . 20 ) 

( 6 . 21 ) 


In equations (6.20) and (6.21) F* is the sonic value of F and and A 
are defined by 
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^i-i/2 

^i+1/2 


F* 


i- 1/2 




F* - F 


,n 


i+ 1/2 


for 

'* 1 - 1/2 

> U* 

for 

“i- 1/2 

< U* 

for 

“i+l /2 

> u* 

for 

“i+ 1/2 

< u* 


(6.22a) 


( 6 . 22 b) 


Both of these schemes produce standard discretizations in regions away 
from sonic lines and supersonic-to-subsonlc shock waves. At sonic lines and 
shock waves these schemes differ from standard schemes and are designed to 
produce smooth solutions through sonic lines and sharp, monotonic shock waves. 
The only difference between the Godunov scheme [eqs. (6.19), (6.20), (6.22)] 
and the Engquist-Osher scheme [eqs. (6.19), (6.21), (6.22)] is in the shock- 
point operator. 


6.4 Free-Stream Consistency Conditions 

Grid-generated irregularities, such as mapping singularities, rapid 
stretching, cell skewness, or grid coarseness, manifest themselves in many 
realistic configurations. Examples of rapid stretching and cell skeimess can 
be found in the grids about multielement airfoils, wlng/fuselage configura- 
tions, wing/fuselage/pylon/nacelle configurations, or even simple airfoils 
when the " 0 " mesh mapping topology is used. In addition, the solution near 
the outer computational boundary almost always consists of extremely coarse 
regions of the grid. Ideally, a stable flow-solver algorithm which can handle 
all of the above mentioned irregularities, yet provide uniform accuracy over 
the entire mesh, is desired. 

The accurate capture of free-stream flow near the outer computational 
boundary where the mesh is quite coarse can be a difficult problem. In a 
formulation that is mapped to the computational domain [see eq. (5.3)], it can 
be shown that if the metric differencing is implemented properly, the trunca- 
tion error associated with a free-stream distribution of the dependent 
variable is zero. That is, free stream is admitted as a solution to the 
finite-difference equations. This type of procedure was addressed by Pulliam 
and Steger (ref. 104) for the Euler equations but was not used, because of the 
small improvements in accuracy obtained on smooth meshes. Thomas and Lombard 
(ref. 105) and Hindman (ref. 106) also worked with geometrically Induced 
errors associated with the metric differencing and found that certain differ- 
encing procedures are better than others. 

All of the above work was performed on the Euler equation formulation. 
Chattot et al. (ref. 107) developed a spatial differencing scheme which con- 
tained a perfect free-stream capture characteristic for the full-potential 
equation. However, in this formulation the full-potential equation was not 
written in strong conservation-law form, that is, the metrics were written 
outside the main flux differentiation as follows (see ref. 108): 

Ai(p(^^)^ + A 2 [(p(^^)^ + (P<t>^)gl + A 3 (p(J>^) + G = 0 (6.23) 

where A^, Ag , and A 3 are defined by equation (5.5) and G is defined by 

G » P-f^ + V^n (6.24) 
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On smooth meshes , where the metric variation is small , this formulation 
behaves like conservative form. On nonsmooth meshes, such as one generated 
by a- solution-adaptive procedure, the metric variation at (for example) a 
shock wave, would be large, this couM greatly affect the conservation prop- 
erties of the spatial differencing stVheme. 


Flores et al. (ref. 109) presented a free-stream-preserving, spatial- 
differencing scheme for the conservative full-potential equation written in 
general curvilinear coordinates [see eq. (5.3)]. Unlike the Euler equation 
scheme presented in reference 104, which produces perfect free-stream capture 
with a single free-stream consistency condition, the full-potential equation, 
in general, requires three conditions. The first condition is associated with 
the density calculation procedure and is developed as follows: The density 

can be written solely as a function of the fluid speed. Thus, the exact 
numerical prediction of free-stream density must result in the exact predic- 
tion of the free-stream fluid speed. The fluid speed can be written as 

^ (6.25) 

This expression reduces precisely to q^, if the difference operators used 
for all ?-dlf ferences involving x, y, and (|) are the same, and if the dif- 
ference operators used for all q-differenceS involving x, y, and (p are the 
same. This can easily be verified by substituting difference operators for 
all derivatives into equation (6.25) and then using the exact free-stream 4> 
distribution to sitapllfy. 

The second and third free-stream consistency conditions are associated 
with the flux calculation. Using the fact that the density is exactly a con- 
stant in free-stream flow, the full-potential equation can be rearranged to 
give 


( ?(h + 5(ti\ /r) Ip +q(fc\ 

^x'*'x ^ ‘^y‘*'yj ^ ^ j „ q 

where and <Py are given by 

, ‘t’x = ^x^ = “ 

‘*’y = S'*"? ^ ^ ^ 


(6.26) 


(6.27a) 

(6.27b) 


If the difference operators for the ^-differences of x, y, and <(> and the 
difference operators for the q-dif ferences of x, y, and iji are the same, 
respectively, then = u„ and i}iy = v„. This can be verified from equa- 
tions (6.27) and is the second free-stream consistency condition. Note that 
this condition is the same as the first condition, providing the density and 
flux calculations are performed at the same grid locations. However, the 
density and flux calculations, in an optimal calculation procedure, are not 
computed in the same locations (see sec. 6.2), and, therefore, these two 
conditions have' to be satisfied independently. 


With 


(jijj = u„ and (|)y = v„, equation (6.26) can be rewritten as 


■ *Sn> ■ " 


(6.28) 


In general, equation (6.28) can be rewritten as two separate equations given 
by 
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- y 


5n 




- X 


en 


- 0 


(6.29) 


Thus, the finite-difference operators used to obtain the metric quantities 
must commute. This Is the third free-stream consistency condition. The last 
condition Is the same condition stated In reference 104 and was required (by 
Itself) to achieve perfect free-stream capture for the Euler equations. 

The three free-stream conditions Just presented are satisfied by the 
two-dimensional finite-volume scheme discussed in section 6.1. However, the 
extension of this finite-volume scheme to three dimensions does not satisfy 
all three free-stream flow consistency conditions (see ref. 76). 


As stated in reference 109, extension of the first two consistency con- 
ditions to three dimensions is straightforward. Extension of the third con- 
dition to three dimensions Is somewhat more difficult. The following three 
equations must be satisfied: 

'V? ' 

- »n"c>{ ^ ‘Vc ■ 

* 'Ve ■ Ve^ * ■ Ve^ • “ 


That is, once all the derivatives of equations (6.30) are replaced by differ- 
ence operators, these relations must cancel Just as in the analytical case 
(see ref. 110) . 


Values of Aj , Aj , A3, and J in the two-dimensional flux calculations, 
satisfying both the second and third consistency conditions, were obtained 
with the same set of difference operators for all primitive metrics, Xjr, x,^, 
etc. (see sec. 5.1). This is not the case in three dimensions. Separate 
formulas for the primitive metrics are required to satisfy the latter two 
conditions. For example, the Aj metric quantity can be written 


A. - CC C ^ + CC 


(6.31) 


The quantities without bars must be computed so as to satisfy free-stream 
consistency condition two, and the barred quantities must satisfy condition 
three. 


The use of this scheme requires more storage (or a moderate Increase in 
execution time if the metrics arc recomputed each iteration). Separate 
values of the metrics are required at four locations, although not all 
metrics are required at each of these locations. The minimum number of 
arrays required for a achem^i in which the density is computed at 
1 + 1/2, J + l/2,k + 1/2 (assuming no metric recomputation) is 15, which com- 
pares to 7 arrays for a standard scheme using the same assumptions. If den- 
sities are computed at i + 1/2, J,k the storage for a perfect free-stream 
capture algorithm Is reduced to 13 arrays. With the size of computer memor- 
ies rapidly increasing, such storage requirements may not be too difficult to 
obtain. 


Two results showing some of the advantages of schemes that capture free 
stream perfectly are shown in figures 19 and 20 (taken from ref. 109). 
Results from a transonic airfoil calculation computed with the TAIR computer 
code (see ref. Ill) are presented in figure 19 for the Kom airfoil at its 
design condition (M» ■ 0.75, a ■ 0.115"; see ref. 112). Three solutions are 
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presented In figure 19 corresponding to (1) a scheme that satisfies none of 
the consistency conditions but uses fourth-order-accutate metric differencing 
(scheme 1), (2) a scheme that satisfies only the first consistency condition 
(scheme 2), and (3) a scheme that satisfies all three consistency conditions 
(scheme 3) . The TAIR code uses a density and metric numerical smoothing at 
the trailing edge of airfoils. For this study, the density smoothing Is 
disabled for all three schemes, and the metric smoothing is disabled for 
schemes 2 and 3. The metric smoothing is still active for scheme 1, for it 
is required for stability. 

In figure 19, the shock location is about the same for all three results. 
However, the shock strength is different in each case. For this calculation 
a shock-free solution is expected, and, therefore, the strength of the shock 
produced is a qualitative estimate of the numerical error associated with 
each scheme. As expected, the third scheme produces the smallest error based, 
on this criterion. Another benefit associated with the third scheme is the 
smooth tralling-edge solution. Because the density smoothing has been 
removed at the trailing edge, schemes 1 and 2 produce oscillations in the 
trailing-edge pressure distribution. However, the scheme 3 result is 
oscillation-free, even without the tralling-edge smoothing terms. Thus, 
satisfying free-stream consistency produces a reduction in global error, as 
well as a reduction in local error around various grid singularities. 

Results from a mesh refinement study are presented in figure 20. The 
lift coefficient is plotted versus the average mesh spacing for the three 
schemes outlined above, all applied to an NACA 0012 airfoil at ■ 0.75 
and a «= 2°. During this study, as the mesh was refined, the ratio of grid 
points on the airfoil surface to the total number of field points was held 
fixed. All three curves, representing the three different types of metric 
differencing, approach the same asymptotic limit as they must to be mathemati- 
cally consistent. On coarse meshes, the error associated with each scheme is 
quite different: scheme 3 is the most accurate, scheme I is next, and 

scheme 2 is the least accurate. The scheme 3 improvement in accuracy on 
coarse meshes represents a highly desirable quality in three-dimensional 
problems. The behavior between schemes 1 and 2 is somewhat unexpected, since 
scheme 2 satisfies the first free-stream consistency condition and scheme 1 
satisfies none of the consistency conditions. Apparently, the fourth-order- 
accurate metrics associated with scheme 1 produce smaller levels of error in 
lift (because the present mesh is smooth) than the second-order-accurate 
metrics of scheme 2. 


6.5 Nonisentropic Full-Potential Formulation 

An interesting exposition of various types of potential formulations 
available for approximating the Euler equations is presented by Klopfer and 
Nixon (ref. 113). Besides the standard mass-energy formulation in which the 
momentum is not conserved, other formulations that conserve momentum and 
energy or mass and momentum are discussed. Crocco's theorem is rederived 
with suitable conservation errors included for each of these formulations. 

It is shown that the Isentroplc assumption is not necessary in conjunction 
with the velocity potential formulation. 

An interesting nonisentropic full-potential formulation is derived in 
reference 113 such that the resultant shock polar is Identical to the Euler 
shock polar. The nonisentropic potential formulation is given by (two- 
dimensional Cartesian coordinates) 

(P^)^ + (P^^)„ - 0 (6.32) 

X A y y 
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p = jK-^ 1 - (^2 + 4>J)]j (6.33) 

where 4> is the standard velocity potential, p is the fluid density, y is 
the ratio of specific heats, and K is a function of the entropy (s) 
defined by 


K = e(s-So)/cv (6.34) 

In equation (6.34) Sq Is the free-stream value of the entropy and Cy is 
the specific heat at constant volume. 

Assuming that the shock waves are normal to the u-component of velocity, 
a locally one-ditnensional assumption can be made (although as pointed out in 
ref. 113, this need not be the case). From equation (6.33) and the Euler 
shock-jump relations, the following relation for K can be derived: 


K 


1 - 


Y 

Y 


Y - 1 

Y+1 


- 1 
+ 1 



(6.35) 


where is the normal velocity upstream of the shock. Values of K equal 

to 1 correspond to free-stream flow; values of K larger than 1 correspond 
to an entropy increase associated with a shock wave. 

Implementation of this scheme requires the ability to detect the posi- 
tion of a shock wave, which is easily accomplished for normal shocks. Then 
the upstream velocity (uj) must be extracted from the flow field. Next, K 
is computed from equation (6.35) for each streamline. This can be done 
iteratively, but in reference 113 a "C" mesh topology is used and it is 
assumed that each ? coordinate line is approximately a streamline. The 
final modification involves a change in the Kutta condition. The far-field 
circulation is no longer given by the velocity potential jump at the trailing 
edge. Now it must be computed from the airfoil circulation and the circula- 
tion generated around the wake, which is due to the jump in K across the 
wake . 


Computed results taken from reference 113 are shown in figure 21 for an 
NACA 0012 airfoil at a free-stream Mach number of 0.80 and an angle of attack 
of 1.25°. The three curves shown correspond to an Euler solution (ref. 33), 
a full-potential solution (ref. 34), and a solution from the Klopf er-Nixon 
nonisentropic potential formulation. For this case, the full-potential solu- 
tion exhibits a strong shock at the airfoil trailing edge. The local shock 
Mach number is approximately 1.5, far exceeding the full-potential formula- 
tion limitation. The nonisentropic formulation, however, essentially pro- 
duces the Euler solution. Even the lower-surface shock wave, which was not 
predicted in the isentropic full-potential formulation, is accurately pre- 
dicted by the nonisentropic formulation. Thus, Euler-like solutions can be 
produced with the nonisentropic full-potential formulation in just a fraction 
of the computer time required by the Euler equations. 

Nonunique solutions to the full-potential equation were reported by 
Steinhoff and Jameson (ref. 114). That is, depending on initial conditions, 
several drastically different flow-field solutions were obtained for the same 
airfoil coordinates with the same free-stream conditions. These nonunique 
solutions only occurred for a range of free-stream Mach numbers that produced 
relatively strong shock waves. Thus, the cause for the nonuniqueness may be 
associated with the large disagreement that exists between the isentropic 




full-potential and Euler shock polars for strong shock waves. If this Is the 
case, the nonlsentropic full-potential formulation of Klopfer and Nixon, 
which possesses the Euler shock polar, may represent a solution to the 
nonunique full-potential problem. 


6.6 Other Spatial Differencing Schemes 


Many other spatial differencing schemes suitable for the solution of 
transonic flow problems based on potential or potential-like foirmulatlonB 
have been presented. A few of these schemes are briefly discussed in this 
section. Of particular note are the field panel method of Piers and Slooff 
(ref. 115); the finite-element method of Vlgneron et al. (ref. 116), and the 
penalty function method described by Bristeau et al. (ref. 117), Periaux 
(ref. 118), and Bristeau et al. (ref. 119). In the last approach, a least- 
square finite-element formulation is used to discretize the full-potential 
equation in conservative form. To exclude expansion shocks the least-square 
functional is modified to include a penalty function. This penalty function 
takes on large values for solutions containing nonphysical expansion shocks, 
that is, for solutions with streamwise positive jumps in velocity, and small 
values for solutions with proper entropy increasing shocks. In a sense, this 
penalty function approach is a dissipative device similar to artificial vis- 
cosity that is designed to exclude physically unrealistic expansion shocks. 


Another approach suitable for the prediction of transonic flow fields is 
based on the stream function formulation. Early pioneering work was done in 
this area by Emmons (refs. 120-122). More recently, Hafez and Lovell 
(ref. 123) presented a stream-function formulation suitable for solving 
transonic flow, x.ie stream-function equation is given by (two-dimensional 
Cartesian coordinates) 



(6.36) 


(6.37) 


and w is the vorticity which is computed from the entropy rise across the 
shock (As). 


In this approach the stream-function equation is discretized using the 
artificial density concept of reference 84 to stabilize the scheme in super- 
sonic regions. A unique solution to the doubled-valued density problem, 
typically associated with transonic applications of the stream-function for- 
mulation [see eq . (6.37)], is presented by Hafez and Lovell. Since the 
irrotationality assumption does not have to be made iJi the stream-function 
formulation, results with vorticity and entropy increases across shock waves 
can be computed. This is highly attractive since solutions to the Euler 
equations can be simulated with this formulation while only requiring the 
numerical solution of a single second-order PDE. In addition, because the 
form of the stream-function equation and the full-potential equation are 
nearly the same, most solution schemes for the potential formulation can be 
used to solve the stream-function formulation. 


The two drawbacks associated with the stream-function formulation are its 
increased complication caused by the double-valued density relation and 


42 


ORIGINAL PAGE IS 
OF POOR QUALfTV 


difficulty in extending It to three-dimensional flows. A three-dimensional 
formulation is presented In. reference 123 but not implemented. The f emula- 
tion expands to three equctlons which are probably still easier to solve than 
the Euler equations; however, until comparisons are made, final conclusions 
cannot be drawn. 


7. ITERATION SCHEMES 


The next subject of discussion involves the iterative process by which 
the initial solution is evolved into the final solution. The iteration 
scheme is primarily responsible for the amount of computational work associ- 
ated with each algorithm through the number of iterations required for con- 
vergence. In recent years, many researchers have experimented with different 
techniques for reducing the number of iterations associated with transonic 
flow computations. In this section we will look at the iteration schemes that 
have proved to be most successful in reducing the computational cost relative 
to older more standard Iterative schemes such as SLOR. 


7.1 The Alternating Direction Implicit Scheme 

The alternating direction implicit (ADI) scheme has already been dis- 
cussed (sec. 3.7) for solving Laplace's equation. In this section, an ADI 
scheme suitable for solving the conservative full-potential equation for tran- 
sonic flow is presented. The ADI factorization used here is basically the 
same as that discussed in section 3.7 and can be stated by writing the stan- 
dard N-operator as follows (see also ref. 108); 


K.i - W.Ki 


(7.1) 


where a is an acceleration parameter (to be discussed shortly) and Aj^ and 
Aj are defined by 


A 


i 



.j 





(7.2) 


In equation (7.2) the density coefficients, p and p are defined by equa- 
tion (6.11) and A^ , A 3 , and J are metric quantities defined by equa- 
tions (5.5) and (5.6). The ADI scheme of equation (7.1) is implemented in a 
two-step format given by the following. 


Step 1: 

(a - ^^A^^^)f" J = uuLiijJ J (7.3a) 

Step 2: 

In equations (7.3), u is a standard relaxation factor and fj is an inter- 
mediate result stored over the entire finite-difference mesh. The residual, 
h<^i,j> is defined by equation (6.10). Step 1 consists of a set of tridiagonal 
matrix equations along the ^ direction, and step 2 consists of a set of 
tridiagonal matrix equations along the n direction. The construction of 
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this ADI scheme does not automatically provide the necessary (|igj. to stabi- 
lize supersonic regions. However, this type of term can still be included 
by adding 




and 


+6 V, J6 
n' i.j' n 


(7.4) 


Inside the parentheses of the first and second steps, respectively. The 
double-arrow notation on the operators of expressions (7.4) Indicate that the 
difference direction is always upwind, and the sign Is chosen so as to 
Increase the magnitude of the matrix diagonal coefficient. The contravariant 
velocity component scaling used in expressions (7.4) provides a smooth tran- 
sition from forward to backward difference directions when the flow direction 
changes sign. The and 6^^ coefficients are constants specified by the 
user according to need. 

Stability of the ADI scheme as given by equations (7.1)-(7.3) can be 
Investigated by considering a simplified form given by 


•(a - 6 )(a - 6 )C? . + ao)(6 +6 )<J>” = 0 

xx"' yy' 1,3 '■ xx yy ^i.j 


(7.5) 


This is essentially the scheme presented in section 3 for solving Laplace's 
equation. Other, more complex model equations can be used, but the present 
one will allow essentially the same conclusions with less work. The amplifi- 
cation factor for the algorithm of equation (7.5) is given by 


®ADI 


^ “ U)35 


+ ^2 


where 


(7.6) 


a 


1 


a + 


4 

aAx^ Ay^ 


(cos aAx 


l)(cos bAy - 1) > 0 


(7.7) 


a 


2 


2 

Ax^ 


(cos aAx - 1) (cos bAy - 1) > 0 

Ay' 


(7.8) 


The purely real amplification factor is always less than or equal to 1 , pro- 
viding 0 < ( 1 ) < 2 and a > 0. Because the only condition for stability on 
the a parameter is that it be positive, the ADI scheme is said to have 
unconditional linear stability, as expected for an implicit scheme. 


The ADI amplification factor can be factored into a special form given 
by (for w = 2): 


2(1 - cos aAx) 

2(1 - cos bAy) 

Ct ““ 

Ax^ 

Ot “ ,r-, -T. ™ . „ : d 

Ay^ 

■ 2(1 - cos aAx)! 

. 2(1 - cos bAy)' 

a + 

Ax^ J 

Ct T* 

Ay^ 


(7.9) 


Note that the x and y directions decouple. Thus, the amplification factor 
can be minimized for a particular eigenvalue associated with (for example) 
the x-dlrection by choosing a to satisfy 

(1 - cos aAx) 


_ 

“ Ax=^ 
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The amplification factor for this value of a and for this eigenvalue Is 
actually zero. For a problem with NI grid points in the x-direction, 
corresponding to NI eigenvalues, the solution will converge to zero error 
after NI Iterations. Of course, this Is only true for a linear problem. 
Precise estimation of the eigenvalues for a nonlinear problem Is generally 
not possible. Instead, a repeating sequence of a's is used with each ele- 
ment of the sequence chosen to maintain small values of |g| for a given 
range of eigenvalues (ref. 124). A suitable sequence of a's presented in 
Ballhaus et al. (ref. 125) is given by 

“k “ k - 1, 2, 3. . , M (7.11) 


where M is the number of elements In the sequence. The sequence endpoints 
can be estimated by using equation (7.10). For Instance, the lowest eigen- 
values, corresponding to low-frequency errors, are approximately given by 
a “ I, which yields 


2 

a B (1 - cos aAx) 

U , 2 

Ax 



B 1 

For high-frequency errors aAx - it, which yields 



(7.12) 


(7.13) 


In practice, it is well advised to "optimize" both and by trlal-and- 
erroT numerical experimentation. Of course, this has to be done only once 
for each code, for and do not strongly depend on the characteristics 
of the solution being computed. 


The ADI scheme just presented has been used to compute transonic flow in 
a number of different applications; however, the results will not be presented 
until after the next section. The Iteration scheme presented in the next , 
section, AF2, will then be compared and contrasted with the ADI scheme in 
section 7.3. 


7.2 The AF2 Approximate Factorization Scheme 


The AF2 scheme was first presented in reference 126 for solving the low- 
frequency (unsteady) TSD equation. This algorithm was subsequently applied 
to the solution of the steady TSD equation (ref. 125) and the conservative 
full-potential equation (refs. 83 and 84). The AF2 fully implicit scheme can 
be expressed by choosing the standard N-operator as follows: 

i «■ - - Vi^eK.J 


where, as with the ADI scheme, Aj^ and Aj are defined by 




> J 


and 



.j-l/2 


(7.15) 
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In equation (7.15) the density coefficients, p and p, are defined by equa- 
tion (6.11) and A 2 > and J are defined by equations (5.5) and (5.6). 

The AF2 scheme Is Implemented In a two-step format given by the following. 

Step 1: 

(a - ^ A )f” . = ao)L(j)” . (7.16a) 

n j i»j i-»j 

Step 2> 

■ 'Id 

In equations (7.16), oi Is the usual relaxation factor, a Is a convergence 
acceleration parameter cycled over a sequence of values [see eq. (7.11)], and 
fj^j Is an Intermediate result stored over the entire finite-difference mesh. 
Step 1 consists of a set of bldlagonal matrix equations along the n-dlrectlon, 
and step 2 consists of a set of trldlagonal matrix equations along the 
^-direction. With the AF2 factorization, the q-dlfference approximation is 
split between the two steps. This generates a i)>r)t~type term, which is useful 
to the iteration scheme as time-like dissipation. The split n term also 
places a sweep direction restriction on both steps, namely, in the negative 
n-direction for the first step [eq. (7.16a)] and in the positive n-direction 
for the second step [eq. (7.16b)]. Flow direction imposes no sweep direction 
restrictions on either of the two sweeps. 

The N-operator, as presented in equation (7.14) (see ref. 34), is some- 
what different from the AF2 scheme presented in references 83 and 125. For 
the AF2 factorization, the N-operator must be written so that either the 
C- or the n-difference operator is split between the two factors. This con- 
struction generates either a (t>gt-type or a 4>^j.-type term and, if properly 
differenced, provides time-dependent dissipation to the convergence process. 
For the "0" mesh topology (see fig. 10), an algorithm with the ^-direction 
split produces a term which is properly differenced either above or 

below the airfoil. Since the supersonic zone can generally exist on both the 
upper and the lower airfoil surfaces and since the supersonic zone is usually 
the most difficult in which to maintain computational stability, it is 
desirable to keep the term differenced in the upwind direction on both 

the upper and lower surfaces. Thus, the N-operator presented in equa- 
tion (7.14) from reference 34 splits the n-direction. This allows control 
of the term because it is added explicitly and is not part of the 

factorization. 

Of course, with this N-operator construction, the term is upwind 

differenced in the forward half of the mesh and downwind differenced in the 
aft half. For this formulation, adverse effects for cases with supersonic 
flow at the trailing edge may be anticipated but none have been experienced. 

In fact, cases with free-stream Mach numbers near un'ty have been computed, 
in which the trailing edge is entirely embedded in supersonic flow with no 
adverse effects (ref. 86). The precise reason for this behavior is unclear. 
The spatial differencing scheme, which is differenced in the upwind direction, 
and the explicitly added term, which is also always differenced in the 

proper manner, may overshadow any adverse effects introduced from the 
term. 

The 4>5t"*^yP® term is included (if necessary) by adding 

+a86^ (7 . 17) 
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inside the parentheses of step 2 [see eq. (7.16b)]. As with the ADI scheme, 
the double-arrow notation indicates that the difference direction is always 
upwind, and the sign is chosen so as to increase the magnitude of the matrix 
diagonal coefficient. The parameter 6 is fixed at a value of 0.3 in sub- 
sonic regions. In supersonic regions, 6 is initialized according to user 
specifications, for example, to 4.5, and then updated using special logic 
(see ref. Ill) given by 


where 


If 

RATIO <2.0 

then 

3 " 

= 0.98 

If 

RATIO >2.1 

then 

3 " 

= 1.1 3 ' 

If 

3 " > BHIGH 

then 

3 " 

= BHIGH 

If 

3 " < BLOW 

then 

3 " 

= BLOW 


RATIO 


ravg" rmax" 
ravg""^ rmax''“^ 


BHIGH = 6^ + 1 , BLOW =6^-1 


(7.18) 


(7.19) 

(7.20) 


In equation (7.19), M is the number of elements in the a sequence, and 
RAVG is the nth iteration average residual. The logic defined by equa- 
tions (7. l8)-(7. 20) monitors solution convergence through the parameter RATIO. 
If convergence is progressing satisfactorily, 6 is reduced; if not, 3 is 
increased. The parameters BHIGH and BLOW are upper and lower bounds which 
limit the amount of 3 variation. This update scheme for 3 is similar to 
the scheme presented by South et al. (ref. 87). 


In addition to the above logic, other larger increases (or decreases) in 
3 are possible. During the iteration process, the base value of 3i includ- 
ing the values of BHIGH and BLOW, are increased or decreased, if the devel- 
oping solution requires more or less time-like dissipation. This logic, 
largely developed on a trial-and-error basis, automatically keys on, the 
growth rates of the number of supersonic points and the amount of circulation. 
If these quantities grow rapidly, then 3, BHIGH, and BLOW are all increased; 
if they grow slowly, then these quantities are decreased. Thus, with this 
type of logic, the time-like dissipation can be automatically adapted to each 
individual solution. 


Stability of the AF2 iteration scheme [eqs. (7 . 14)-(7 . 16) ] can be inves- 
tigated by considering a simplified algorithm given by 


(at - 6 )(a - ■? )C" + au(6 + 6 )(]." . 

' X yy''' x' i,j xx yy i,J 


= 0 


(7.21) 


This is essentially a scheme for solving Laplace's equation. The standard 
von Neumann test yields an amplification factor given by 


where 


a^ + (1 - bi)a^ 

°AF2 “ a^ + a2 


(7.22) 


a, = ^ (1 - + ? (cos biy - l)(e^^^’' - 1) (7.23) 

* aAxAy^ 
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— ^ (cos bAy - 1) i 0 (7.2A) 

Ay^ 


Unlike the ADI scheme amplification factor, which was purely real, the AF2 
amplification factor is complex and, therefore, a little more difficult to 
analyze. However, the imaginary parts in both the numerator and denominator 
are the same. Thus, an equivalent condition of stability is given by 


-I < 


aj + (1 - w)a2 


a, + a- 


< 1 


(7.25) 


where 


aj = Re(a^) 

= (1 ~ cos aAx) + (cos bAy - l)(cos aAx - 1) > 0 (7.26) 

oAxAy^ 

The expression given by equation (7.25) is of the same form as the amplifica- 
tion factor for the ADI scheme [eq. (7.6)1 and, therefore, is always satis- 
fied, providing 0 < u < 2 and a > 0. Thus, the AF2 scheme has the uncon- 
ditional linear stability typically associated with a fully implicit 
iteration scheme. 

There is an interesting aspect of the stability of the present AF2 
scheme associated with the airfoil surface boundary condition, as discovered 
by South (ref. 127). The residual operator airfoil surface boundary condi- 
tion is that of flow tangency and is implemented by reflection. Another 
aspect of the surface boundary condition is that the n-direction difference 
on f at the airfoil surface is arbitrarily set to zero. South noted that a 
stability analysis using the proper airfoil surface boundary conditions (with 
reflection condition included) produces the following stability condition: 

a > pAjO) (7.27) 

That is, the a parameter (or equivalently the inverse of the time-step) is 
restricted at the airfoil surface by the condition given in equation (7.27). 
The quantity A 3 is effectively the cell aspect ratio (see sec. 5.1). Thus, 
as the mesh is clustered tovjard the airfoil surface, the stability condition 
(7.27) becomes more restrictive and, unless the a sequence is suitably modi- 
fied, divergence can result. 

South presented a solution to this problem which consisted of an appro- 
priate modification of the N-operator at the airfoil surface. By modeling 
the surface N-operator after the surface residual operator, unconditional 
linear stability was restored. 

Finding values of a that cause the AF2 amplification factor to be zero 
for a given eigenvalue is difficult because the AF2 amplification factor is 
complex. However, values of a can be determined which do minimize the AF2 
amplification factor (see ref. 125) and are given for the* low- and high- 
frequency limits as follows: 


“l " ’ “h “ (7.28) 

Notice that the high-frequency endpoint for the AF2 scheme is quite different 
from that of the ADI scheme but that the low-frequency values are the same. 
Again, in practice, it is well advised to optimize both endpoints by 
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trial-and-error numerical experimentation. This, however, need only be done 
once for each code, for solution parameters such as Mach number, angle of 
attack, and airfoil coordinate variations do not greatly affect the optimal 
a sequence. 


The AF2 algorithm discussed above has been coded Into a user-oriented 
computer code called TAIR (Transonic AIRfoll analysis). Computational 
results produced with TAIR are now presented. The first result Involves the 
supercritical Kom airfoil at a free-stream Mach number of 0.74 and an angle 
of attack of 0°. Pressure coefficient distributions for this slightly off- 
deslgn case are compared In figure 22, with a result from the GRUMFOIL com- 
puter code (ref. 128). The GRUMFOIL computer code has available a viscous 
correction option which was not used for this calculation. The two results 
are in excellent agreement. The rms error convergence history curves for 
this calculation are presented In figure 23. The rms error at iteration 
n (E^g) Is defined by 


where Cp^ 
the nth 


RMS 


fNI 


E 

(CPi - Cp^)^ 

i=i 



NI 




1/2 


(7.29) 


is the surface-pressure coefficient at the ith grid point and 
iteration; Cp^ is the surface-pressure coefficient, at the ith 


grid point taken from the converged solution; and NI is the total number of 
surface grid points . Using ErjiJs to compare convergence performance is a 
much more quantitatively correct procedure than using the standard maximum 
residual quantity. (More discussion of this point can be found in refs. 83 
and 129 and in sec. 7.3 of these notes.) The three curves shown in figure 23 
correspond to the following iteration schemes: (1) AF2, (2) hybrid, and 

(3) SLOR. The hybrid scheme is a combination semidirect/SLOR iteration 
scheme developed by Jameson (ref. 85), which is composed of one semidirect- 
solver iteration followed by several SLOR iterations. The purpose of the SLOR 
iterations is to smooth high-frequency errors generated by the direct-solver 
step in regions of supersonic flow. For this calculation, the AF2 scheme con- 
vergence rate is about 5 times faster than that of the hybrid scheme and about 
10 times faster than the SLOR rate. 


Additional results obtained from the TAIR computer code are shown in 
figures 24 and 25. Figure 24 shows a pressure coefficient comparison with 
experiment taken from reference 111. The airfoil is the supercritical CAST 7 
and the Mach number is 0.7. The agreement is quite good. Figure 25 shows 
Mach-number contours around an NACA 0012 airfoil immersed in a 0.95 Mach free 
stream at an angle of attack of 4° (taken from Holst, ref. 130). The Mach 
number contours clearly illustrate the existence of a so-called "fishtail" 
shock-wave pattern downstream of the airfoil trailing edge. This difficult 
calculation demonstrates the convergence reliability associated with the AF2 
transonic relaxation procedure. 

The AF2 scheme has been implemented in a number of different applica- 
tions, including airfoil calculations by Holst (ref. 34) and Atta (ref. 131); 
cascade flows by Kwak (ref. 132); wing geometries by Holst (ref. 130) and 
Holst and Thomas (ref. 133); and in wing/pylon/nacelle calculations by Atta 
and Vadyak (ref. 134). In addition, two other factorizations similar to the 
AF2 scheme, have been presented by Benek et al. (ref. 135). These two formu- 
lations are called AFZ2 and AFZ and are both used for three-dimensional 
transonic wing calculations. The AFZ2 scheme is very similar to the standard 
AFZ scheme used in reference 130. The AFZ scheme is a simplification of the 
AFZ2 scheme which inverts matrices along only the wraparound and normal-like 
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directions, not the span direction. Thus, this scheme is Implicit In only 
two directions, whereas the standard AF2 scheme Is Implicit along all three 
coordinate directions. Despite this feature, the simpler AFZ scheme has 
almost the same convergence properties as the AFZ2 scheme. Numerical results 
using the three-dimensional AF2 scheme of reference 133 will be presented in 
section 8. 


7.3 Convergence Characteristics of SLOR, ADI, and 
AF2 Iteration Algorithms 

Numerical results comparing the convergence characteristics of the two 
fully Implicit algorithms just presented (ADI and AF2) with the SLOR algo- 
rithm are now presented. All three iteration schemes have been applied to 
the same artificial-density, spatial-differencing scheme for the conservative 
form of the full-potential equation. A two-dimensional, 10%-thlck, circular- 
arc airfoil with small-disturbance boundary conditions is used as a test case. 
The finite-difference grid is Cartesian with variable spacing in both the 
X and y directions. Both subcritical and supercritical cases are considered 
(M„ = 0.7 and 0.84, respectively). Pressure coefficient distributions for 
these two cases are displayed in figure 26. Note the perfect symmetry asso- 
ciated with the subcritical case and the existence of a moderate strength 
shock at about 80% of chord for the supercritical case. For more details 
about these calculations see reference 83. 

Convergence characteristics for the subcritical case are displayed in 
figure 27. All of the convergence parameters for each scheme have been 
selected by a trial-and-error optimization process. Based on a six-order-of- 
magnltude reduction in the maximum residual, the ADI scheme is about twice as 
fast as the AF2 scheme and about 16 times faster than SLOR, in terms of iter- 
ation count. However, the ADI and AF2 schemes take about 50% and 30% more 
CPU time per iteration, respectively, than SLOR; this should be considered 
when speed ratios based on the total amount of computational work are desired. 

Convergence characteristics for the supercritical case are displayed in 
figure 28. Again, the convergence parameters have been optimized by a trial- 
and-error process. Based on a six-order-of-magnitude reduction in the maximum 
residual and in terms of iteration count, AF2 is slightly more than twice as 
fast as ADI, and about 11 times faster than SLOR. The number of supersonic 
points (NSP) plotted versus iteration number for the supercritical case is 
shown in figure 29. The AF2, ADI, and SLOR schemes reach the final value of 
NSP in 29, 103, and 320 iterations, respectively. 

The AF2 iteration scheme was relatively consistent in terms of conver- 
gence speed for both cases. The ADI iteration scheme, on the other hand, dis- 
played remarkable speed for the subcritical case but was a disappointment for 
the supersonic case. This is because the <l>st“*^yP® error term produced by the 
AF2 factorization is more suitable for supersonic regions than the ((ij^-type 
error term resulting from the ADI factorization. In fact, the (Ji^-type error 
term has been shown to be destabilizing in the supersonic region (see 
ref. 29). 

Relative levels of convergence for AF2 and SLOR for given reductions in 
maximum residual are compared in figure 30. The solid lines represent the 
final solidly converged solution. The other results represent intermediate 
AF2 and SLOR solutions in which the maximum residual has been reduced by one, 
two, and three orders of magnitude (figs. 30(a), 30(b), and 30(c), respec- 
tively) . It is immediately obvious that reducing the maximum residual by 
equal amounts for the AF2 and SLOR schemes does not produce intermediate 
results with the same level of error.' This behavior can also be observed by 
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comparing the maximum residual history curves of figure 28 with the rms error 
history curves given in figure 31. The rms error is computed from a formula 
similar to equation (7.29). The SLOR residual drops very rapidly initially 
and then levels off. The SLOR rms error drops gradually. Therefore, at the 
"knee" in the SLOR residual history curve, even though the residual has 
dropped by about three orders of magnitude, the actual rms error has dropped 
by only one order of magnitude. In contrast, both maximum residual and rms 
error results for the ADI and AF2 schemes are nearly straight lines with 
about the same slope. 

This behavior is the result of two factors (see refs. 83 and 129): 

(1) the AF2 scheme treats all error components equally well (approximately), 
whereas the SLOR scheme performs efficiently on only the high-frequency error 
components; and (2) it can be shown that the residual is a weighted sum of 
errors, in which the weighting factors are the eigenvalues of the finite- 
difference scheme. The eigenvalue for high-frequency errors is 0(Ax~ ); for 
the low-frequency errors it is 0(1). Hence, the residual is heavily influ- 
enced by the high-frequency errors. Therefore, the maximum residual operator 
should not be used as the basis for comparing convergence performance between 
iteration schemes with different characteristics (for example, AF and SLOR 
schemes). The rms error is much better suited for this purpose. In prac- 
tice, using the maximum residual to monitor convergence for either AF or 
SLOR is the most convenient method (since error is unknown). However, the 
convergence criterion based on residual should be adjusted (by experience) in 
accordance with the solution procedure in use. 


7.4 AF3 Approximation-Factorization Scheme 

Another interesting approximate-factorization scheme, introduced by 
Baker (ref. 136) for solving the nonconservative full-potential equation is 
now discussed. The AF3 scheme is given by 

(-aCly - A<5^j^)(a + = au)L4>"^j (7.30) 

for subsonic regions of flow and by 


(-aC t - A 6 
' c y c XX 


^n 


- A ■^ ) (a + t )C" . - au)L<|.V , 
u x-^ ' y' i,j i,j 


(7.31) 


for supersonic regions of flow (when u > 0) . Ihe difference direction on the 
third term of the first factor is reversed when u < 0. The coefficients A, 
C, Ac, Cj., and Ay are determined from the nonconservative full-potential 
equation written in canonical form (stream and stream-normal coordinates) , 


and are given by 



1 


(7.32) 


A = Ay + Aj, , 


C = Cy + Cc J 


where u and v are the velocity components along the x and y directions, 
respectively; a is the local speed of sound; and q is the magnitude of the 
velocity vector. As in the previous AF2 scheme, m is a relaxation factor. 
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a an acceleration parameter (cycled over a sequence o£ valuea), and . 

is the nth iteration residual defined by 


b<t> 


i.J 


(A6 + By u 6 6 + C6 . 

^ XX *^x^y X y yy'^l.j 


for subsonic regions of flow and by 


(7.33) 
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(7.34) 


for supersonic regions of flow (when u > 0 and v > 0) . The additional 
coefficients By and Bj. are given by 



The AF3 factorization is similar to the AF2 factorization and produces 
exceptional improvements in computational speed. This fact is illustrated in 
figures 32 and 33 (taken from ref. 6). Figure 32 shows the pressure coeffi- 
cient distribution about an NACA 0012 airfoil at «• 0.75 and a •• 2“ for 
several different levels of convergence. After just 10 AF iterations, which 
correspond to about 13 SLOR iterations, the solution is nearly converged. 

The AF3 convergence history for this case is compared with an SLOR convergence 
history in figure 33. The SLOR convergence history is enhanced by using the 
standard grid refinement procedure Involving two grids, one coarse and one 
fine. For this case, the SLOR scheme required over 400 iterations to reach 
plottablc accuracy; the AF3 scheme reached plot table accuracy in about 
20 Iterations (26 equivalent SLOR iterations). 


The two-dimensional AF3 scheme just discussed has been extended to 
three dimensions by Baker and Forsey (ref. 137). Solutions of the nonconser- 
vative full-potential equation have been obtained for wing and wing/fuselage 
combinations with a factor of 4 or 5 increase in computational efficiency 
relative to standard SLOR schemes. 


A theoretical analysis for various AF schemes, including ADI, AF2, and 
AF3, is presented by Catherall (ref. 138). In this study, an improvement in 
computational efficiency is obtained when the contributions of the transfot- 
matlon metrics are properly split between the two factors. In addition, 
optimal values for the acceleration parameter sequence a and the relaxation 
factor ui are' derived and discussed. 


7.5 Multigrid Iteration Schemes 

The multigrid iteration scheme is enjoying a wave of popularity that has 
included applications in a host of different areas. This scheme is actually 
a convergence acceleration technique and requires a base iteration scheme, 
for example, SOR, SLOR, or AF. Multigrid schemes have existed for quite some 
time, having been first introduced by Fedorenko (ref. 139) in 1964. Since 
then, several authors have analyzed the technique, including Bakhvalov 
(ref. 140) and, more recently, Nicolaides (ref. 141) and Hackbusch (ref. 142). 
The most significant aspect of the multigrid iteration scheme is fast con- 
vergence. This fast convergence is produced by using a sequence of grids 
ranging from very coarse to very fine. Each grid is used to eliminate one 
small range of errors in the error frequency spectrum, namely the errors of 
highest frequency supported on each mesh. Many relaxation schemes exist that 


work very well on high-frequency errors, for example, point-Jacobi and AF 
schemes (with properly chosen acceleration parameters). One of these relaxa- 
tion schemes is used on each mesh to remove the high-frequency error. A 
desirable aspect of this approach is that the high-frequency error on the 
coarsest mesh is actually the lowest frequency error existing in the problem. 
Because this usually troublesome low-frequency error is efficiently dealt 
with on a coarse mesh, very little computational work is expended in removing 
it from the solution. Thus, a tremendous convergence rate enhancement is 
obtained . 

Implementation of a typical multigrid scheme is described in general 
terms as follows: Suppose we desire a solution to 

L% = f (7.36) 

where is a typical linear difference operator which approximates a dif- 

ferential operator L on a mesh associated with the grid spacing h. The 
quantity f contains the problem boundary conditions. Let 

(J) a. u + V (7 . 37) 

where u is an approximation to <|) and v represents an error. Therefore, 
as the iteration scheme converges, u ->■ and v ->• 0. The basic multigrid 
scheme can be expressed by 

L®^v + - f) = 0 (7.38) 

n 

9 Vi 

where L is a finite-difference operator which approximates L on a mesh 
associated with the grid spacing 2h, instead of h, that is, twice as coarse 
as the original mesh. The operator is an interpolation or averaging 

operator which transfers values of the residual (L u - f) from the fine mesh 
to the coarse mesh. After the coarse mesh corrections, v, are obtained, they 
are transferred back to the fine mesh by using 


= u + I^.v (7.39) 

k 

where l 2 h interpolation operator. The process can continue to 

coarser meshes so that ultimately just one or maybe several mesh cell widths 
span the entire domain of interest. 


To extend the idea to nonlinear problems a simple modification is help- 
ful. By adding and subtracting L^^u from equation (7.38) the new form 
becomes 

= I (7.40) 


where 


u = u + V 

f = L^^u - - f) 

h 


(7.41) 


The quantity u represents a new or improved estimate of (|i which is deter- 
mined from the coarse mesh. The quantity f is a modified right-hand side 
which essentially represents the difference in residuals between the h and 2h 
meshes. New updated coarse values are transferred back to the fine mesh by 
using 
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(7.42) 


u”®” U + 

Thus, the error quantity v does not have to be stored as in the original 
version. 

Applications utilizing the multigrid scheme were slow to materialize 
after its introduction, primarily because of difficulties in Implementation 
and general underestimation of the potential of multigrid enhanced schemes. 

The first work to apply the multigrid scheme numerically was that of Brandt 
(ref. 143) in 1972. Later, the multigrid scheme was formulated in general 
terms by Brandt (ref. 144). In this latter reference, a good historical 
background of the multigrid scheme is presented, including a review of 
related earlier work. 

The first use of the multigrid scheme for transonic calculations was 
presented by South and Brandt (ref. 145). In that study, numerical solutions 
of the TSD equation for nonlifting airfoils were obtained. The speed of an 
optimized SLOR scheme was improved by a factor of 5 on uniform meshes and by 
a factor of 2 on stretched meshes. A primary difficulty reported by South 
and Brandt involved the existence of a variety of limit-cycle oscillations 
between several grids, thus inhibiting convergence. This problem seemed to 
be the result of insufficient smoothing of the high-frequency errors on one 
grid before passing to the next coarser grid. South and Brandt concluded that 
the SLOR base algorithm used in the raultigr5.d scheme did not have uniform 
smoothing properties in both directions, especially for nonuniform, highly 
stretched meshes. They hypothesized that alternating the SLOR sweep direction 
or utilizing an ADI iteration scheme as the base algorithm might solve this 
problem. 

Another approach, proposed by Arlinget (ref. 146), is to refine or 
coarsen the mesh in only one coordinate direction while doing line relaxation 
along the opposite direction. This technique produced a convergence rate 
acceleration but did not take full advantage of the multigrid scheme. To 
date, the most successful application of a multigrid convergence acceleration 
scheme to a practical transonic problem is the work of Jameson (ref. 147). 

In that study, the full-potential equation in conservative form is 
solved, using a multigrid scheme with a specially constructed AF base itera- 
tion scheme. This scheme, when applied to the following linear model 
equation. 


Ad) + Bd) 
^xx ^yy 


(7.43) 


is given by 


(S - 


(7.44) 


where A and B are constants, oi is the standard relaxation fector, and 
S and L are operators defined by 


S “ On -f a,6 +0,6 

o lx 2 y 


(7.45) 


Lii" . (a6 ^ + B6 )(|)" . 

^i,j XX yy ^i,j 
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In equation (7.45) a^, and are parameters that depend on flow type 
and user input. The quantities 6~ and fiy denote first-order-accurate 
upwind difference operators in the x and y directions, respectively. 

The Jameson scheme uses a recursive approach for implementing the multi- 
grid philosophy. Instead of the adaptive approach advocated by Brandt 
(refs. 143 and 144). In the adaptive multigrid approach, the decision to 
proceed to the next mesh, either coarser or finer, is based on a convergence 
rate criterion. If the solution residual is dropping slowly, the Iteration 
process proceeds to coarser meshes. Conversely, if the solution residual is 
dropping rapidly, the iteration proceeds to finer meshes. In the recursive 
approach of Jameson, a single multigrid cycle starts with an AF iteration on 
the finest mesh, followed by an AF iteration on the second finest mesh, etc. 
This continues until the coarsest mesh la reached. Then the process is 
reversed, starting with the coarsest mesh and ending with the second finest 
mesh. Therefore, each multigrid grid cycle consists of one AF application on 
the finest mesh and two applications on each of the remaining meshes. If a 
fine grid AF iteration is defined as a unit of work, then one multigrid 
cycle, using the recursive approach, requires about 1-2/3 work units plus 
interpolation operations. 

Results produced by the Jameson multigrid scheme are displayed in fig- 
ures 34 and 35. The pressure coefficient distribution for an NACA 64A410 air- 
foil at a free-stream Mach number of 0.72 and at an angle of attack of 0“ is 
displayed in figure 34. A moderate-strength shock exists at about 60% of 
chord. Notice that the residual has been reduced below lOE-12 (see fig. 34), 
which is approximately an eight-order-of-magnitude reduction from the initial 
value; the reduction was achieved in only 29 multigrid cycles. Convergence 
histories for this case, which were computed using different numbers of 
meshes (from one mesh, that is, no multigrid, up to five meshes), are shown 
in figure 35. The convergence rate (CR) , which is defined as the mean 
reduction in the average residual per unit of work, is also displayed for 
each curve. Increasing the number of meshes or, equivalently, increasing the 
coarseness of the coarsest mesh, greatly improves the convergence rate. 

Other researchers have used the multigrid algorithm to sove the full- 
potential equation in a variety of applications: Fuchs (ref. 148) and 

Deconinck and Hirsch (ref. 150), for two-dimensional applications; Ar linger 
(ref. 150), for axisymmetric calculations; McCarthy and Rehner (ref. 151) and 
Brown (ref. 152), for three-dimensional engine-inlet calculations; and 
Shmilovich and Caughey (ref. 153) and Caughey (ref. 154), for three- 
dimensional wing calculations . 


7.6 Other Iteration Schemes 

The strongly implicit procedure (SIP) Introduced by Stone (ref. 155) has 
been applied to the numerical solution of the full-potential equation for 
transonic airfoil calculations by Sankar and Tassa (ref. 156). Additional 
applications include those of Sankar et al. (ref. 157), for steady transonic 
wing calculations, and Roach and Sankar (ref. 158), for transonic cascade 
calculations. In all cases the SIP solution algorithm displayed good conver- 
gence characteristics as a relaxation scheme. In addition, the SIP algorithm 
has the ability to compute time-accurate flow fields; see Sankar et al. 

(ref. 159) for unsteady wing calculations. 

The SIP iteration scheme requires three additional arrays of storage 
(five arrays for the SIP scheme, two for the ADI or AF2 schemes) and requires 
a few more operations to invert the resulting matrix equations. The SIP 
method has a built-in mechanism for matrix conditioning. That is, the scheme 
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automatically adjusts the matrix diagonal entries so that the matrix is well 
conditioned, even in regions where the mesh spacing becomes coarse. Addi- 
tional properties regarding the SIP iteration scheme are discussed in 
reference 157. 

Several researchers have presented new algorithms especially designed 
for vector computers, including Keller and Jameson (ref. 160), Hafez et al. 
(ref. 84), Redhead et al. (ref. 161), Hotovy and Dickson (ref. 162), and 
South et al. (ref. 87). Vector computers offer greatly enhanced computing 
speeds arising from the ability to operate on many calculations simultaneously 
(parallel machines) or in an assembly line fashion (pipeline machines). 

[See Bailey (ref. 163) for more discussion of computer architectures.] In 
references 84, 160, and 162 the algorithm's vector characteristics were 
stressed above all else. As a result, these vector algorithms could process 
very large numbers of grid points per second but usually required more 
iterations than standard SLOR to converge. 

In South et al. (ref. 87) an algorithm called Zebra II, which is highly 
vectorj.zable and requires about the same number of iterations to converge as 
does SLOR, is described. This algorithm is an explicit, or point, scheme 
which mimicks a full-plane SOR algorithm and is designed to solve the conser- 
vative full-potential equation. 

The Zebra II algorithm takes a step in the right direction, but other 
approaches for vector computation may still be superior. A theoretical study 
comparing the vector processing attributes of four transonic full-potential 
algorithms (SLOR, ADI, AF2, and Zebra II) was performed by Holst (ref. 164) 
utilizing a mathematical model of a pipeline vector computer. The results of 
this study indicate that implicit algorithms, which contain nonvectorizing 
matrix inversions, still enjoy an overall supremacy on vector computers rela- 
tive to explicit, or point, iterative techniques, when all aspects of effi- 
ciency are taken into consideration. Results from a vectorization study per- 
formed on a three-dimensional transonic wing code will be presented in the 
next section. 

Other iteration schemes suitable for producing fast convergence for the 
full-potential equation include the extrapolation schemes of Hafez and Cheng 
(ref. 165), Caughey and Jameson (ref. 49), and Yu and Rubbert (ref. 166); the 
conjugate-gradient methods of Bristeau et al. (ref. 119), Glowinski et al. 

(ref. 167), Chattot and Coulombeix (ref. 168), and Wong and Hafez (refs. 169 
and 170); and the minimum residual method of Wong and Hafez (ref. 171). The 
work of Wong and Hafez (ref. 169) provides an interesting discussion of itera- 
tion schemes for solving the full-potential equation. Results are presented 
for several schemes, including SLOR, two variations of the ZEBRA scheme men- 
tioned above, and conjugate-gradient schemes with several types of precondi- 
tioning combined with both SLOR and ZEBRA schemes. These schemes are applied 
to two spatial discretization schemes, including a finite-difference scheme 
and a finite-element scheme. It is found that the combined iteration schemes 
are superior to the standard SLOR scheme by as much as a factor of 10 for 
subcritical cases and by at least a factor of 2 for tough transonic cases. 
Details of each of these schemes can be found in the references cited above. 

The multiline, or block, iterative schemes presented in Hafez and Lovell 
(ref. 172) represent another competitive form of relaxation algorithm. In 
this type of scheme two or three lines of the grid are treated in the same 
matrix inversion. Thus, the amount of implicitness normally associated with a 
single line inversion scheme such as SLOR, is greatly increased. Several 
convergence rate comparisons are presented in reference 172 and indicate com- 
petitive computational efficiencies for the multiline schemes relative to 
other types of iteration schemes. In addition, with red-black ordering of the 
blocks, these schemes are easily vectorized. 
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8. APPLICATIONS IN THREE DIMENSIONS 


Extension of the ideas presented in the previous chapters to three 
dimensions, for example, wing and wlng/fuselage calculations, is very much 
complicated by storage and computer time requirements. Nevertheless, the 
benefits to be obtained are large and much work has been done in this area. 

In this section a brief review of existing techniques including specific 
aspects of a number of widely used three-dimensional transonic computer codes 
is presented. For a more detailed review of this and other related subjects, 
see reference 2, 


8.1 Computer Code Characteristics 

Several computer codes designed to solve the three-dimensional transonic 
flow over wing and vjing/fuselage configurations have been developed and are 
currently in wide use in the aircraft industry. A few of these codes are 
listed in table 2 along with specific references describing the details of 
each code. Note that both conservative and nonconservative forms of the 
full-potential formulation are represented. More discussion concerning the 
various transonic potential formulations, including both full-potential and 
TSD forms, is given in chapter 2 of reference 2. 

The iteration schemes used by the codes listed in table 2 are predomi- 
nantly SLOR. The newer approximate-factorization or multigrid schemes are 
utilized in only the more recent codes. This produces very slow and thus 
costly convergence for most of the older codes and makes them relatively 
expensive to run. One feature generally used in SLOR codes to improve com- 
putational efficiency is the use of grid refinement, that is, the use of a 
coarse-medium-fine grid sequence to accelerate convergence. Converged results 
from the coarse mesh are interpolated onto the medium mesh, then reconverged 
and finally interpolated from the medium mesh to the final fine mesh. Thus, 
a good initial solution is provided for the fine-mesh calculation. Use of 
this grid sequence philosophy increases the computational efficiency by a 
factor of 2 or 3 (at least for crude levels of convergence). 

The overall computer time required for a complete transonic wing solution 
varies greatly, depending on the level of convergence desired, the number of 
grid points used, ind the type of formulation chosen. Even the amount of 
computer time reported by two different users from the same computer code can 
be different by factors of 3-5, simply because of the level of convergence 
desired and the number of grid points used. Computer time comparisons which 
address most of these aspects are given in table 3 (taken from Hinson and 
Burdges (ref. 185). Note that all codes compared have about the same number 
of wing-surface grid points, whereas the total number of field points is quite 
different. The Bailey-Ballhaus code (see ref. 186) uses a concept called grid 
embedding (first introduced by Boppe (ref. 187) to mote adequately cluster 
grid points at the wing surface) . This concept is very attractive because 
the resulting mesh topology is much more efficient and the wing-surface 
results do not seem to be affected by this treatment. Grid embedding has not 
been used extensively with the full-potential formulation, probably because 
of the resulting mapping complications. In addition, full-potential grids 
are generally more efficient than TSD grids (that is, TSD grids without 
grid embedding), and therefore, grid embedding would not be as useful in the 
full-potential formulation. (See Brown, ref. 188, for one application using 
the full-potential equation and embedded grids.) 

All three results shown in table 3 were obtained for the same transonic 
wing at identical test conditions. All three codes used a coarse-medium-fine 
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sequence. Convergence was monitored by examining the pressure distribution 
histories at two span locations. The Bailey-Ballhaus TSD code converged in 
approximately 200 iterations, the FL022 code in about 50-100 iterations, and 
the FL027 code in about 200 iterations. For these conditions, the FL027 code 
solution time is approximately twice that of the FL022 code and about 5 times 
that of the Bailey-Ballhaus code. The convergence characteristics were 
essentially the same for both the conservative and nonconservative options 
available in the Bailey-Ballhaus code. For more information about this set 
of comparisons see reference 185. 

Another computer time comparison is reproduced in table 4 from refer- 
ence 183. In this comparison, results from the FL028 computer code on the 
CDC 7600 computer and the TWING computer code on both the CDC 7600 computer 
and the new Cray-lS vector computer are presented. The recently introduced 
TWING computer code is similar to the FL028 code in that it solves the three- 
dimensional full-potential equation in conservative form, but different in 
that it uses the fully implicit AF2 iteration scheme (see sec. 7.2). 

The results for both codes presented in table 4 have been obtained from 
the same problem (ONERA M6 wing, = 0.84, a = 3.06°). The convergence 
parameters from both codes have been approximately optimized by a trial-and- 
error process. The TWING results are based on an 83-iteration run in which 
the center-span lift changed by 0.03% in the last 30 iterations. The FL028 
results are based on a coarse-medium-fine mesh sequence with 50, 50, and 
283 iterations, respectively. The total lift for this run changed by 0.49% in 
the last 100 iterations of the fine mesh. Because FL028 uses a less efficient 
mapping, more field points are required to achieve approximately the same num- 
ber of surface grid points used in the TWING code. 

Immediately obvious from table 4 is that TWING sets up the lift much 
faster than FL028 (14 times faster for 98% of the lift) when both codes are 
run on the same computer (CDC 7600). Of course, this improvement is due to 
the advanced fully implicit iteration scheme utilized by TWING (see sec. 7.2 
for more discussion of this point). In addition, the vectorized version of 
TWING run on the Cray-lS vector computer is about 11 times faster than TWING 
run on the CDC 7600 computer. These run times are for the flow-solver portion 
of both codes only and do not include times for grid generation, initializa- 
tion, or solution printout. 

The speed improvement offered by TWING is, of course, largely a result 
of the faster computer hardware associated with the Cray-lS computer. But 
another important reason for the improvement is that the AF2 iteration scheme 
in the TWING code is highly vectorizable. All the Interior grid point opera- 
tions performed by the TWING code vectorize, including all operations asso- 
ciated with the matrix inversions of each sweep, the upwinding logic asso- 
ciated with the density coefficients, and the logic associated with the 
time-like damping term. Despite this complete vectorizatlon, the code is 
almost completely written in FORTRAN and will run on any standard computer 
with a FORTRAN compiler. The FL028 computer program would benefit from the 
vectorizatlon offered by the Cray computer, but would have difficulty obtain- 
ing the level of efficiency achieved by the TWING computer code. This is 
because a portion of the SLOR algorithm used within FL028 is inherently recur- 
sive, and, therefore, not vectorizable. For more Information regarding these 
results see reference 183. 

The Jameson-Caughey nonconservative full-potential computer code (FL022) 
uses a sheared parabolic coordinate system which can be defined by (see 
ref. 173) 
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= tx “ Xq(z) + l[y - y^(z)J}^^^ 

Zi = z 

where z is the spanwlae coordinate, and Xq and y^, define a singular line 
of the coordinate system located just Inside the leading edge. The effect of 
this transformation is to unwrap the wing to form a shallow bump (see fig. 36): 


( 8 . 1 ) 


Yi = s(xi,yi) 

Next, a shearing transformation is used 

C * X 

n = y^ - s(xi,yi) 

C = z 


( 8 . 2 ) 


(8.3) 


to map the wing surface to a coordinate surface. Finally, the h, and C 
coordinates are stretched with suitable stretching formulas. A more reveal- 
ing diagram of the overall wing grid system shown in the physical domain is 
reproduced from reference 189 in figure 37. Note that the airfoil cross- 
sectional grid produced by this transformation technique is of the "C" mesh 
variety. 


An interesting aspect of the FL022 algorithm, created by the grid topol- 
ogy, is the orientation of the tridiagonal matrices used in the SLOR iteration 
scheme. This aspect is illustrated in figure 38, which has been taken from 
reference 173. This diagram shov;s the computational domain in a spanwise 
cross-sectional plane; note the orientation of the streamlines. The purpose 
of sweeping through the mesh in this fashion is to always avoid sweeping in an 
upwind direction in the supersonic region. If this procedure is not followed, 
unstable operation could result. 


The James on-Caughey FL027 computer code is capable of treating the flow 
about isolated swept wings or wings mounted on an infinite cylinder. The 
FL028 code is a closely related derivative of FL027 which allows more sophis- 
ticated treatment of the fuselage. The FL030 code allows a still more sophis- 
ticated treatment of the wing/fuselage problem. The relative level of 
sophistication of each of these formulations is displayed in figure 39 (taken 
from ref. 190). The actual transformation used for most FL027 and FL028 wing 
calculations is very similar to the transformations described for the FL022 
computer program and will not be discussed further. 

The ONERA transonic full-potential code of Chattot et al. (ref. 107) is 
capable of treating the flow about isolated wing geometries. The mapping 
used is similar to the previous descriptions and will not be discussed further. 
This code, like the more recent FLO codes (FL027, FL028, and FL030) , uses 
numerical evaluations of the metric quantities and, therefore, could support 
mesh generation routines of varying types provided they all used the same 
general topology. 


Most full-potential codes are designed for either conservative or noncon- 
servative spatial differencing schemes. However, the ONERA code includes 
options for both forms. This allows direct comparison of the two forms in a 
format that is unencumbered by different iteration algorithms, grid topologies, 
or programming styles. Such a comparison is shown in figure 40 for a rectan- 
gular planform, nonlifting NACA 0012 wing at a free-stream Mach number of 
0.85. The conservative result produces a shock that is too strong and about 
10%-15% downstream of the experimental shock position. The nonconservative 
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shock location is incorrect in just the opposite direction; it is too weak 
and too far forward. This type of disagreement between conservative and non- 
conservative forms is characteristic. The code user should be aware of which 
formulation is being used and how the resulting solutions should be 
interpreted. 

The ARA full-potential wing/body code of Forsey and Car (ref. 178), is 
capable of solving the nonconservative full-potential equation about isolated 
wings and wing/body combinations. Bodies must be circular in cross-section 
although the body radius may vary axially. This has the effect of providing 
a relatively good wind/body interference simulation but does not produce good 
finite-length body effects. 

The mapping procedure used in the ARA code is implemented in two phases. 
First, the body is transformed via conformal mapping into a slit that corre- 
sponds to a portion of the symmetry plane. Because each body cross section 
is constrained to be circular, this mapping is analytic. Next, the wing lead- 
ing and trailing edges are extended to infinity to produce a "flat-plate wing" 
outboard of the tip. Each flat-plate wing cross section is mapped to the 
interior of a circle using an analytic conformal mapping, and each finite- 
thickness wing cross section is mapped using a numerical conformal mapping. 
Thus, the total effect of this transformation is to map the surface of the 
wing (including the flat-plate wing extension) to the inside surface of a 
circular cylinder and the free-strearo outer boundary to the circular cylinder 
center. Unlike the FLO full-potential codes (F1022, FL027 , FL028, and FL030) , 
which use a "C" mesh topology about each wing span station, the ARA code uses 
an "0" mesh topology. This produces a somewhat more efficient mesh for the 
ARA code; that is, for the same number of total grid points, the ARA code mesh 
puts about 30% to 40% more points on the wing surface than a typical FLO code 
mesh. 


The TWING code is capable of treating isolated-wing or simplified wing/ 
fuselage geometries. The finite-difference grid is generated by an elliptic- 
solver numerical grid-generation procedure; it is described in Holst and 
Thomas (ref. 133). An example grid generated using this procedure is dis- 
played in figure 41 (taken from ref. 191). The geometry used is a highly 
swept, low-aspect-ratio wing with a large amount of twist (Hinson-Burdges 
Wing C, see ref. 185). Figure 41(a) shows the entire planform view of the 
wing. Three airfoil cross-section plots of the grid are displayed in fig- 
ures 41(b)-41(d): (1) at the wing root, (2) just inboard of the tip, and 

(3) just outboard of the tip, respectively. The exact positions of these 
cross-sectional plots are shown in the planform view of figure 41(a). The 
x/c and z/c Cartesian coordinates displayed in figures 41(b)-41(d) are 
plotted to the same scale and are normalized by the root chord, c. Note the 
large amount of twist and the efficiency with which the "0" mesh topology 
clusters grid points about each airfoil cross section. 

Finally, figure 41(e) shows a perspective view of the Wing C grid as 
generated by the TWING computer code. The symmetry plane, wing surface, and 
vortex sheet grids are all displayed. For clarity, only every fourth grid 
line in the wraparound direction is plotted. This view (except for the 
trailing-edge vortex sheet) corresponds very closely to the view of the actual 
wing mounted in the wind tunnel. 


8.2 Results Obtained with Existing Codes 

The next item of importance is accuracy. That is, how well do these 
codes actually predict transonic-wing, flow-field physics? In this section. 



computed results collected from a variety of sources will be displayed in an 
attempt to quantify the answer to this question. 

A series of results taken from reference 107 is shown in figure 42. 

Three numerical results, Jameson (conservative), Bailey /Ballhaus (TSD con- 
servative, see ref. 186), and ONERA (nonconservative) are compared with 
experiment at a wing semispan station of 20%. The configuration used for 
this set of calculations is the ONEPA M6 wing at M„ ■= 0.841, a = 3°. (See 
ref. 192 for more information on the experimental results and geometrical 
characteristics of the ONERA M6 wing.) This frequently used ONERA M6 test 
case represents a difficult test for a transonic flow calculation procedure 
because of the existence of a double shock. The double-shock structure is 
preducted by all three results with some variation in the second shock posi- 
tion. The conservative results predict a somewhat stronger second shock 
slightly downstream of the nonconservative shock. 

Results from the ARA full-potential code (taken from ref. 2) are pre- 
sented in figures 43 and 44. The wing/body configuration used in this cal- 
culation is shown in figure 43. The body was a circular cylinder of constant 
radius, and the mid-mounted wing had the following characteristics: AR = 6, 

TR = 1, sweep = 25°. The wing had a constant supercritical airfoil section 
and was used without twist. Viscous effects were simulated by the addition 
of a displacement thickness obtained from a two-dimensional transonic viscous 
code (see ref, 193). 

Computed results for a free-stream Mach number of 0.86 and an angle of 
attack of 4.2° are compared, in figure 44, with experimental results for the 
same Mach number and two angles of attack, 4.6° and 3.7°. The calculation 
was performed with the usual grid sequence involving 300, 100, and 160 itera- 
tions on the coarse, medium, and fine meshes, respectively. The finest mesh 
for this case consisted of 160 x 20 x 24 = 76,800 grid points. The computa- 
tional time required for this calculation was equivalent to about 25 min of 
CPU time on the CDC 7600 computer. Results are shown for three span stations 
(y = 0.37, 0.55, and 0.73), and are generally in excellent agreement with 
experiment. Note the unusually sharp shock capture displayed in figures 44(b) 
and 44(c). This is a result of the efficient mesh topology used in the ARA 
full-potential code in conjunction with a relatively fine mesh. 

Additional results computed from the ARA full-potential computer code 
are displayed in figures 45-47 (taken from ref. 137). Figure 45 shows the 
wing/fuselage geometry used in this calculation. For more details regarding 
this configuration, which has been designated as an AGARD test case for three- 
dimensional method evaluation, see reference 194. Pressure distributions for 
this geometry at Mo, = 0.9 and a = 1° are displayed in figure 46. These 
results were computed on a relatively coarse mesh consisting of 
80 X 24 X 10 = 19,200 points. The agreement between the computed and experi- 
mental results is good at all three semispan stations shown in figur.'i 46. 

The results just presented have been computed with the new version of the 
ARA code which uses the AF3 approximate-factorization scheme of Baker and 
Forsey (ref. 137), Lift development versus the equivalent number of relaxa- 
tion iterations is presented for this case in figure 47. Four curves are 
displayed corresponding to AF and SLOR iteration schemes applied to the wing/ 
body configuration of figure 45 and a wing-alone variation of the same con- 
figuration. The AF3 convergence rate displayed in figure 47 is far superior 
to the SLOR convergence rate for both the wing-alone and the wing/fuselage 
cases. The SLOR scheme used in figure 47 did not use the standard grid 
sequence which would have improved convergence. However, the AF3 scheme still 
possess a large advantage over the SLOR scheme in convergence speed. 



Pressure comparisons for the FL022 and FL027 computer codes are pre- 
sented in figures 48 and 49 for a high-aspect-ratio (AR = 10.3) supercritical 
wing tested by NASA Langley (taken from ref. 195). The wing is swept 27° at 
the quarter chord line, and has streamwise section thickness ratios of 14.9% 
at the fuselage junction, 12% at the traillng-edge break, and 10.6% at the 
tip. (See ref. 196 for a more detailed description of the geometry and 
experimental results.) Viscous effects were modeled by a two-dimensional 
integral boundary-layer method applied in streamwise strips. The calculation 
angles of attack were determined by matching the experimental lift 
coefficients. 

Figure 48 presents results from FL022; FL022 does not model the fuselage 
used in the experiment . In an attempt to account for the fuselage interfer- 
ence effect, this calculation was run at a free-stream Mach number of 0.80, 
whereas the experiment was run at 0.79. Agreement is good everywhere except 
near the inboard stations, where details of the fuselage interference are 
probably important. 

Figure 49 presents results from FL027 for the same case as that shown in 
figure 48. For this case, however, the fuselage is modeled as an infinite 
cylinder with a radius equal to the maximum radius of the fuselage on which 
the wing was mounted in the experiment. The low-mounted position of the test 
wing was also simulated in the numerical calculation. For this case, the 
aeroelastic deformation of the model was estimated by introducing 0.36° of 
wash out (negative twist) at the tip. This aeroelastic twist was added 
linearly from the trailing-edge break to the tip; it was not included in the 
FL022 calculation of figure 48. The Mach-number correction used for the 
FL027 result was 0.007 and was applied to account for the effect of the 
finite fuselage used in the experiment. It is clear from figure 49 that 
representation of the effect of a finite fuselage by a simple Mach-number 
shift is an oversimplification. The resulto in figure 49 are in better agree- 
ment with experiment near the root station than the results of figure 48, but 
in poorer agreement outboard near the tip. Nevertheless, reasonable agreement 
is obtained in both cases by making these ad hoc adjustements or calibrations. 

A relatively simple, but adequate way in which to introduce fuselage 
effects into isolated-wing codes is discussed in Verhoff and O'Neil (ref. 190) 
and by Van der Vooren et al. (ref. 197). In the latter approach, spanwise 
velocities are prescribed at the vertical wing-root plane in the FL022 com- 
puter code. These velocity components are determined from a panel-method 
calculation for the complete wing-body configuration. Pressure distributions 
showing the success of this technique are presented in figure 50. Numerically 
computed results with and without this root-plane modification are compared 
with experiment at four semispan stations. In this calculation the fuselage 
effects are extremely important immediately at the root plane but die out 
rapidly as midspan is approached. 

The next result presented in this section is from the Yu computer code 
(refs. 181 and 182). This code is basically a modified version of the 
Jameson-Caughey finite-volume code. The primary improvement in the Yu code 
is a new efficient and flexible grid-generation method based on the body- 
fitted coordinate method of Thompson et al. (refs. 36-38). In this approach, 
a set of nonlinear elliptic equations is solved numerically to determine the 
three-dimensional grid. The Interior grid distributions are controlled by the 
boundary grid distributions in a manner developed by Thomas and Moddlecoff 
(ref. 39). For more information about this grid-generation technique see 
section 5.2. 

Results from the Yu code (taken from ref. 2) are presented in figures 51 
and 52. The wing/fuselage configuration is shown in figure 51. In this 
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calculation procedure the restriction to simplified fuselage geometries is 
not required; reasonably complex geometries can be accurately simulated. In 
this case, the flow field about the Boeing 747-200 wlng/fuselage geometry is 
computed for a free-stream Mach number of 0.84 and an angle of attack of 2.8”. 
The section pressure coefficient distributions are compared with experiment 
at three span stations In figure 52. The agreement at each station Is 
excellent . 

A comparison of FL028 and TWING results with experiment for the 
ONERA M6 wing at = 0.84 and a = 3.06“ is presented in figures 53-56 
(taken from ref. 133). The finite-difference grid used by TWING for this 
calculation consisted of 89 x 25 x 18 = 40,050 points (wraparound, spanwise, 
and normal-like directions, respectively) with 89 x 17 = 1513 points on the 
wing surface. The FL028 results were computed on a mesh with 
120 X 16 X 28 = 53,760 points (wraparound, normal-like and spanwise direc- 
tions, respectively) with 75 x ig = 1350 points on the wing surface. The 
TWING cross-sectional grid uses the "0" mesh topology, and FL028 uses the 
"C" mesh topology. This produces a somewhat less efficient mesh for FL028 
relative to TWING and is the basic reason for comparing solutions computed 
on different size meshes. 

Pressure coefficient comparisons are presented in figure 53 at four dif- 
ferent semispan stations (n - 0.20, 0.44, 0.65, and 0.90). The general 
agreement of the results is good, especially considering the coarseness of 
the meshes involved. The largest discrepancy arises at the aft shock. The 
FL028 shock is in better agreement with experiment, being very close or 
slightly upstream, and the TWING shock is slightly downstream. Since these 
two codes are conservative inviscld formulations it is expected that the 
inclusion of viscous effects would move this shock forward. Such a movement 
would bring the TWING result into better agreement with experiment and have 
the opposite effect for the FL028 result. 

The overall shock sonic-line position for all three results is compared 
in figure 54. These results correspond to the shock sonic line computed by 
linear interpolation and do not necessarily reflect the position of the 
steepest shock pressure gradient. The experimental points almost exactly 
split the two numerical results , with FL028 being upstream and TWING 
downstream. 

The convergence properties of the TWING computer code for the ONERA M6 
wing problem Just presented are displayed in figures 55 and 56. The conver- 
gence parameters for this case have been approximately optimized by a trial- 
and-error process. Figure 55 shows the buildup of the number of supersonic 
points (NSP) and the center-span, section-lift coefficient with iteration 
number (n) . Also plotted along the horizontal axis is the CPU time asso- 
ciated with both the CDC 7600 and Cray IS computers. (Note that the computer 
times for this case have already been presented and compared with the FL028 
computer times in section 8.1; see table 4). Each symbol plotted (fig. 55) 
represents one cycle in the a acceleration parameter sequence. At 
24 iterations, or three a cycles, the solution is nearly converged. This 
corresponds to just 4.8 sec of CPU time on the Cray IS vector computer! 

The rate of development of the pressure coefficient distribution at one 
selected semispan station (n = 0.18) is presented in figure 56, The solution 
a*'ter 20 iterations is compared with the tightly converged solution after 
80 iterations. Very little disagreement exists, except at the shock, where 
about six or eight points have not quite converged to their final steady- 
state values. 
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The next result presented is a TWING calculation (taken from ref. 191) 
for the Wing C configuration described in Hinson and Burdges (ref. 185) (see 
also ref. 198). This wing geometry is characterized by a small aspect ratio 
(AR “ 2.6), large sweep (Al^ = 45°), small taper ratio (TR ■ 0.3), and has a 
large amount of twist (twist =8°). The grid for this configuration has 
already been presented in figure 41; it consists of 127 x 27 x 20 ■* 68,580 
points with 127 x 17 => 2159 points on the wing surface. Pressure distribu- 
tion comparisons between the TWING results and experiment are shown in 
figure 57. For this case, the experimental Mach number and angle of attack 
were 0.85 and 5.9°. The computational Mach number and angle of attack were 
0.83 and 5.0°. The numerically computed corrections (AM = -0.02 and 
Aa = -0.9°) were determined by trlal-and-error; they represent (approxi- 
mately) the best experimental/numerical porrelation for this set of condi- 
tions. The set of corrections determined for these data in reference 198 
(AM = -0.005 and Aa = -0.9°) are similar to but smaller than those computed 
by TWING. 

The agreement for this case is quite good everywhere except at the tip, 
where the need for viscous corrections is apparent. Of particular note in 
this calculation is the ability of the TWING code to predict the oblique 
shock, which exists at both the third (y/c = 0.5) and fourth (y/c * 0.7) 
semispan stations. The differencing scheme in this region is entirely first- 
order accurate and yet little shock smearing is exhibited. 

To obtain an idea of the effects of the angle of attack (Aa) and the 
Mach number (AM) corrections, several results are shown in figure 58 (taken 
from Subramanlan et al., ref. 191). Pressure distributions at two semispan 
stations for the Wing C configuration just discussed are compared with 
experiment for several conditions; (1) the uncorrected experimental condi- 
tions (M„ =0.85, a = 5.9°); (2) the corrected experimental conditions, 
using the corrections cited in reference 198 (M^o = 0.845, a = 5.0°); and 
(3) the corrected experimental conditions, using the corrections computed in 
reference 191 (M<„ = 0.83, a = 5.0°). As seen in figure 58, the angle-of- 
attack correction is more important than the Mach-number correction. The 
corrections cited in reference 198 yield a reasonable solution in the present 
case, primarily because both angle-of-attack corrections are the same. 

The last set of results is from a version of the FL030 computer code in 
which a sophisticated viscous correction model has been added by Streett 
(ref. 199). This model consists of a three-dimensional. Integral boundary- 
layer method based on the work of Stock (ref. 200) and Smith (ref. 201). In 
addition, a strip wake model accounting for thickness and curvature effects 
is also Included. Results using this computational technique are displayed 
in figures 59 and 60. • 

Figure 59 shows a pressure coefficient distribution compared with 
experiment for the Hinson-Burdges Wing A geometry at Moo = 0.819 and 
a = 1.96° (see ref. 185). Three semispan stations are shown, p = 0.28, 0.50, 
and 0.68. In this case, the experimental Mach number and angle of attack 
are matched. Note the excellent agreement with experiment. 

Figure 60 shows the lift reduction owing to viscous effects, R, which 
is defined by 

c — c 

p = &,inviscid i. 

^£,inviscid 

plotted versus span station for the advanced transport configuration pre- 
sented in reference 196. The Mach number and angle of attack for both the 
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experiment and the computation were 0.78 and 1.65°, respectively. The two 
curves displayed in figure 60 correspond to the full wake model of Streett 
(ref. 199) and a similar wake model without curvature effects. Viscous 
effects are very important for this calculation, as shown by both curves in 
figure 60. In addition, the wake-curvature effect is shown to be substantial 
outboard of the mid semispan position, especially near the wing tip. 

The preceding results are offered in an attempt to answer the question 
posed at the beginning of this section: How well do these schemes model 
transonic flow-field physics? In some cases, the modeling is quite good, in 
other cases not so good. The subjects of geometrical modeling and viscous 
corrections play the most important roles in obtaining good experimental/ 
numerical correlations. Other effects, such as wind-tunnel-wall interfer- 
ence, tunnel flow angularity, model distortion, numerical truncation error 
(coarse mesh effects), and lack of solution convergence, also must be con- 
sidered in the final analysis. 

Most of the three-dimensional results presented in this section show at 
least one new or different attribute; for example, computational efficiency, 
geometrical generality, or good viscous corrections. In most codes, only 
one area is fully developed; thus, many codes exist that are, for example, 
very accurate but very slow or that are very fast but do not contain accurate 
viscous effects. Consequently, there will have to be additional research in 
consolidation of techniques before truly useful transonic codes can be made 
available. The ultimate transonic code must have good characteristics in all 
respects — accuracy, speed, reliability, and the capability of solving flow 
fields over a variety of different aircraft configurations. 


9. SUMMARY AND RECOMMENDATIONS FOR FUTURE RESEARCH 


The numerical .solution of the transonic full-potential equation has 
received much attention within the CFD research community in the past 
10-12 years. The purpose of these notes has been to review selected topics 
associated with this field of research. Classical relaxation schemes and 
early ideas associated with transonic potential schemes were touched upon 
first. Current research acticities were then discussed, with emphasis on 
grid-generation techniques, spatial differencing schemes, and iteration 
schemes. Special emphasis was placed bn the convergence acceleration char- 
acteristics of iteration schemes since this aspect largely controls compu- 
tational cost. Finally, a series of three-dimensional results from a variety 
of different sources was presented, thus providing the reader with a good 
basis for evaluating the state of the art in this field. 

One obvious conclusion from this presentation is that the numerical solu- 
tion of the full-potential equation is highly developed. The aircraft 
industry utilizes, on a routine basis, numerical solutions of the full- 
potential formulation to gain insight into the transonic flow fields for 
many two- and three-dimensional applications. It is projected that with the 
cost of computations going down and the cost of wind-tunnel tests going up, 
the use of these computational tools for the design and development of modem 
aircraft will increase steadily in the future. 

Before the application of these schemes can become truly production- 
oriented, there will have to be improvements in several areas. First, com- 
puter codes must be more flexible and, therefore, applicable to a larger 
range of geometries. Advances in this area are largely paced by developments 
in surface-geometry representation and grid generation. Second, accuracy 
must be improved. This includes relatively simple areas including 
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boundary-condition application, metric differencing, and grid-point cluster- 
ing in regions of large solution gradients, as well as more complicated 
aspects associated with viscous effects. Third, in the area of viscous 
effects, the accurate prediction of moderate separation regions must be 
Improved. And fourth, the computational efficiency and reliability of these 
schemes should be Improved. That is, application of a specific technique to 
a slightly different geometry or new set of free-stream conditions should 
not result in divergence or a significant loss of speed. 

As indicated in the last section, significant progress has been made in 
many of these areas. However, most of the computer codes in existence today 
have serious failings in one or more of the areas mentioned above. Therefore, 
a clear-cut line of future research is in consolidating the advances that 
already exist. The ultimate computer code must contain all of the character- 
istics mentioned in the preceding paragraph: generality, accuracy, reliabil- 

ity, and computational efficiency. 

This type of research is more atuned to applications rather than basic 
algorithm development, but research associated with basic algorithm develop- 
ment'" is also required. Many of these topics are associated with grid- 
generation algorithms. Interface differencing procedures associated with 
overlapping grid systems or component adaptive grids need to be investigated. 
This approach seems to be the beat way (at least for finite-difference dis- 
cretization schemes) of eventually producing transonic flow solutions about 
complicated geometries — for example, over a complete aircraft. Other 
research opportunities associated with the full-potential equation exist in 
the area of solution-adaptive grids. This would produce improvements in both 
computational and storage requirements and at the same time be directly 
applicable to other equation formulations. 

Additional basic research associated with the full-potential equation in 
the areas of computational efficiency or algorithm accuracy or both will 
still be important in the years to come, but the trend in this regard is seen 
to be more in the direction of the Euler and Navier-Stokes equations . With 
these equation sets, more complete and, therefore, more accurate descriptions 
of the flow field physics, will be possible. The main advantage this pro- 
vides is the analysis of off-design conditions with significant viscous 
effects, including massive separation (Navier-Stokes equations) or in the 
analysis of flow fields with significant vorticity Interaction, for example, 
in a closely coupled wing-canard configuration (Euler equations). The major 
problems associated with these equation sets are obtaining sufficient compu- 
tational efficiency, geometrical generality, and, in the case of the Navier- 
Stokes equations, adequate turbulence models. 

In this regard the full-potential formulation may provide a significant 
benefit. Formulations that effectively simulate the Euler formulation (for 
example, the nonisentropic full-potential formulation presented in sec. 6.5) 
can be developed as a much cheaper alternative to the Euler equations. In 
another type of application, the standard full-potential formulation could be 
used to accelerate convergence for the Euler or Navier-Stokes formulations. 

In many transonic cases of interest, solutions of the full-potential and 
Euler formulations are very similar. A full-potential solution could be used 
as an "initial guess" for the Euler calculation, and thus significantly 
reduce the computational effort associated with the more expensive Euler cal- 
culation. In still another application, the full-potential equation could be 
used in regions of the flow field where the irrotationality and isentropic 
assumptions are valid, leaving the regions with significant entropy and vor- 
ticity gradients for the Euler equations. For many calculations these 
regions may be very important but confined to only a small portion of the 
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entire flow field. Thus, a large Improvement In computational efficiency 
can be gained by the u/3e of such a hybrid scheme. 

In summary, research opportunities associated with the full-potential 
formulation are numerous and varied. They range from the development of 
application codes to basic algorithm research. The use of these methods in 
practical applications to predict transonic flow fields about general three- 
dimensional configurations Is Increasing and Is expected to help aircraft 
designers produce more efficient aircraft at less cost. 
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TABLE 1. CONVERGENCE RATE ESTIMATES FOR 
VARIOUS RELAXATION SCHEMES (FROM 
REF. 9) 


Algorithm 

Number of iterations^ 

Point-Jacobi 

2/A^ 

Point -Gauss-Seidel 

l/A^ 

SOR 

1/2A 

Line-Jacobi 

1/A^ 

Line-Gauss-Seidel 

1/2A^ 

SLOR 

1/2 *^A 

ADI 

-log(A/2)/1.55 


“dumber of iterations required for a one- 
order-of-magnitude reduction in error 
[Num, see eq. (3.45)]. 
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TABLE 2. SOME EXISTING THREE-DIMENSKMJAL FULL-POTENTIAL CODES 


Code name 
(ref.) 

bate 

Eq^ 

type“ 

Iteration 

scheme" 

Grid 

Remarks (ref.) 

FL022 

(29,173) 

1974 

NC 

SLOR 

Sheared 

conformal 

Isolated wings, 
viscous option (174) 

Reyhner 

(175,151) 

1976 

NC 

SLOR or 
MG/SLOR 

Sheared 

Cartesian, 

nonaligned 

Axisymmetrlc Inlets 
(n > 0) with or 
without center bodies 

FL027,FL028, 

FL030 

(52-54) 

1977 

C 

SLOR 

Sheared 

conformal 

Wing and wlng/body 
[tunnel walls (176)] 
[AF version (135)] 

Dassault 

(177) 

1977 

C 

Optimal 

control 

conjugate 

gradient 

Finite 

element 

Wing, wlng/body and 
wlng/body /nacelle 

ONERA 

(107) 

1978 

C or 
NC 

SLOR or 
AF 

Sheared 

conformal 

Isolated wings 

ARA 

(178) 

1978 

NC 

SLOR 

Sheared 

conformal 

Wing and wing /body 
[AF version (137)] 

Eberle 

(179) 

1978 

C 

SLOR 

Finite 

element 

Isolated wings 

Chen and 
Caughey 
(180) 

1979 

C or 
NC 

SLOR 

Sheared 

conformal 

Axisymmetrlc inlets 
(a > 0) with and 
without centerbodles 

Yu 

(181,182) 

1980 

C 

SLOR 

Elliptic 

solver 

Wing/body and 
wing/body /nacelle 

Sankar, Malone, 
Tassa (157) 

1981 

C 

AF 

Sheared 

conformal 

Isolated wings 
[Unsteady option (159)] 

TWING 

(133,183) 

1982 

C 

AF 

Elliptic 

solver 

Wing and wlng/body 

TUNA 

(184) 

1982 

C 

AF 

Elliptic 

solver 

Isolated wings 
(Unsteady option) 


* nonconservative; C » conservative. 


SLOR ■ successive-line over relaxation; MG “ multigrid; AF ■ approxi- 
mate factorization. 


TABLE 3. GRID FEATURES AND COMPUTING TIMES FOR SEVERAL TRANSONIC WING 

COMPUTER CODES (FROM REF. 185) 


Code 

Wing chordwlse 
grid points 

Wing spanulse 
grid points 

Total 

field 

points 

Approximate 
CPU time,** 
min 

Number of 
iterations 

Balley- 

Ballhaus 

37 

25 

41,000 

6 

200 

FL022 

61 (root) 
40 (tip) 

21 

159,000 

15 

50-100 

FL027 

51 

21 

90,000 

30 

200 


^CDC 7600 computer. 


TABLE 4. A COMPARISON OF CODE EXECUTION SPEED FOR THE 
ONERA M6 WING AT M„ - 0.84, a ■ 3.06° 

(FROM REF. 183) 



FL028 

TWING; 
scalar 
(ref. 133) 

TWING; 
vector 
(ref. 183) 

TWING: 
vector 
(ref. 183) 

Computer 

CDC 7600 

CDC 7600 

CDC 7600 

Cray- IS 

Field points 

53,760 

40,050 

40,050 

40,050 

Surface points 

1,350 

1,513 

1,513 

1,513 

Time for 98% 

742 

64 

53 

4.8 

lift, sec 
Ratio 

155 

13 

11 

1 



M<1 


M<1 



M>1 

Fig. 1. Typical transonic flan 
field c^out an airfoil. 



Fig. 2. Typical transonic flow 
field about a swept wing. 




Fig. 4. Characteristics and the TSD 
equation. 


FULL POTENTIAL FLOW 



Fig. 5. A comparison of the isen- 
triprio shock-dump relation and the 
Euler shock-dump relations (ref. 12). 


Fig. S. The domains of dependence and 
influence for steady supersonic flow. 
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SIMPLE UPWIND SUPERSONIC DIFFERENCING 
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SLOR 
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ADI 



O DIFFERENCING MOLECULE USED FOR L 
^ IRESIDUAL OPERATOR) 

4- DIFFERENCING MOLECULE USED FOR N OPERATOR 


CHARACTERISTICS 



Fig. 8. Miamatoh between the phyeioal 
and nmerioal domccine of dependence — 
a deatdbiUzing aituation. 


Fig. 6. A summary of various 
olassiaal relaxation schemes. 
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Fig. 7. A comparison of convergence 
rate estimates for the relaxation 
schemes of figure 6. 


Fig. 9, A schematic representation 
of the rotated differencing scheme of 
Jameson (ref. 29). 
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X 


(a) Phys-iaal domain. 


VORTEX 

WEET 

(LOWER) 



(b) Computational domain. 

Fig. 10. Nimerioally generated air- 
foil transformation: (x,y) ->• 

"0" mesh topology. 



Fig. 11. Finite-difference grid about 
a highly cambered 12 to 1 ellipse^ no 
control (ref. 41). (a) Entire geome- 

try. (b) Trailing-edge close-up. 



Fig. 12. Finite-difference grid about 
a highly cambered 12 to 1 ellipsei 
with control terns activated (ref. 41). 
(a) Entire geometry, (h) Trailing- 
edge close-up. 


til 


m 



(a) Entire geometiy. 



(b) Close-up of the trailing edge. 


Fig. 13. Finite-difference grid about 
a highly cambered thin ellipse gener- 
ated via the parabolic scheme of 
Edhamura (ref. 44). 
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Fig. 14. A wing / fuselage grid generated via the parabolia scheme of 
Nakamura (ref. 46). 



Fig. 16. An example solution-adaptive grid applied to a transonic airfoil 
calculation: NACA 0012 airfoil, = 0.75, a = 2° (ref. 76). 
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Fig. 17. Pressure coefficient distributions for the MCA 0012 airfoil at 
= 0.75 and a = 2°, both density and artificial viscosity parameter \> 
computed at node points (ref. 86). 



Fig. 18. Pressure coefficient distributions for the NACA 0012 airfoil at 
= 0.85 and a = 0°, artificial viscosity parameter v computed at node 
points (ref. 87). 
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V SCHEME 3 
•: ;j SCHEME 2 
« O SCHEME 1 

1.0 - 


■ 1,21 1 1 1 1 

•70 .80 .90 1.00 1.10 

x/c 


Cloa^-up of tvailing-edge solution. 


(a) Entire airfoil. (b) 

Fig. 19. Pressure coefficient distributions for the Kom airfoil at its 
design point CM„ = 0.75 and a = 0.12°) showing the effects of different 
spatial differencing schemes on the solution (ref. 109). 



Fig. 20. Mesh refinement study for the NACA 0012 airfoil at M^, = 0.75, 
a = 2°, showing the effects of grid coarseness on several spatial differenc- 
ing schemes (ref. 109).' 
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EULER (STEGER) 

FULL POTENTIAL (HOLST) 
STRONG SHOCK FULL POTENTIAL 
(KLOPFERNIXON) 




TAIR 

QRUMFOII. 


F^g. 21. Avvfo%l pressure aoeffiaient 

comparisons of the nonisentropio and Fig. 22. Pressure coefficient oompar 

isentropia full-potential formulations ison for the Korn airfoil at 

and the Euler formulation (ref. 113). = 0.74 and a - 0° (ref. 111). 


O AF2 
□ HYBRID 
a SLOR 



80 120 
CPU TIME, MC 


Fig. 23. Convergence history comparison for the Kom airfoil at = 0.74 
and a = 0° (ref. 111). 






x/c 


Fig. 26. Suhoritiaal and superaritiaal solutions used for the aonvergenoe 
history comparisons of figures 27-31 (ref. 83). 





Fig. 27, Maximum vesidual aonvergenae history aomparieon, subaritioal case, 
= 0.7 (ref. 83). 



Fig. 28. Maximum residual convergence history comparison, supercritical case, 
= 0.84 (ref. mh 



Fig. 29. Development of the supersonic region, W* - 0.84 (ref. 83). 
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(a) Maximum residual reduced one order of magnitude. 



(b) Maximum residual reduced two orders of magnitude. 

Fig. 30. Intermediate solution comparisons after specified reductions in the 
maximum residual^ supercritical case^ Af„ = 0.84 (ref. 83). 
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ORiGlNAL 

OF POOR QUALITY 



(a) Maximum residual reduced three orders of magnitude. 
Fig. SO. Concluded. 



Fig. SI. Root-mean-square error convergence history comparison, supercritical 
case, = 0.84 (ref. 8SJ . 
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Fig. 32. Devolopment of the pi'eseure aooffioient distvibution with iteration 
nwiibei' /bj’ the /If 3 itemtion sahme^, MCA 0022 au'foil^ » 0.?S, a = 2" 

( ref. B) . 
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ORIGINAL PAG£ ^ 
OF POOR QUALITY 


NACA 0012 
M.-0.7B, a- 2° 



0 100 200 300 400 


EQUIVALENT RELAXATION ITERATIONS 

Fig. 33. Convergence history compari- 
son showing the AF3 and SLOR iteration 
schemes, NACA 0012 airfoil^ = 0.75. 
a = 2° (ref. 6). 



NACA 64A410 

MACH 0.720 ALPHA 0.000 

CL 0.6667 CD 0.0030 CM -0.1478 

GRID 192X32 NCYC 29 RES 0.323E-12 

Fig. 34. Converged pressure coeffi- 
cient distribution from the FL036 
muttigrid code, NACA 64A010 airfoil, 

= 0.72, a = 0° (ref. 147). 



WORK 

NACA64A410 

MACH 0.720 ALPHA 0.000 
RESID1 0.720E04 

GRID 192X32 

Fig. 35. Multigrid convergence his- 
tories, NACA 64A410 airfoil, = 0.72, 
a = 0° (ref. 147). 


t. 



Fig. 36. Construction of the sheared 
parabolic coordinate system used in 
FL022 (ref. 173). 
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Fig. 37. Sheared parabolia coordinate 
system used in FL022 (ref. 189). 


STREAMLINE 


VERTICAL LINE 
RELAXATION 



Fig. 40. Pressure coefficient com- 
parisons for a rectangular plan form 
NACA 0012 wing, conservative versus 
nonconservative differencing 
(ref. 107). 


HORIZONT/^L LINE 
RELAXATION 


MARCHING 

DIRECTION 


Fig. 38. Marching directions of 
relaxation scheme used in FL022 
(ref. 173). 


FLO-27 GEOMETRY 


FLO-28 GEOMETRY 


FLO-30 GEOMETRY 


Fig. 39. Comparison of the various 
FLO code geometry modeling capabili- 
ties (ref, 190). 


(a) Wing plan form. 

Fig. 41. Numerically generated finite 
difference mesh about the Hinson- 
Burdges Wing C configuration 
(ref. 191). 
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(b) Ving-root aross-sectional grid 
(127 X 20 grid points). 


(d) Cross~seotionat grid just out- 
board of tip (127^ 20 grid points). 



(a) Cross -seationaZ grid just inboard 
of the wing tip (12? >• 20 grid points). 



(e) Perspective view of the grid, 

127 X 27 X 20 grid points (only every 
fourth point plotted in the wraparound 
direction) . 


Fig. 41. Concluded. 
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OF 'i*OOR 


JAMESON (C) v/b'0.20 



Fig. 42. Invisaid pressure aoeffi- 
aient comparisons for the OFERA M6 
wing configuration^ = 0.841^ 
a = 3“ (ref. 107). 



Fig. 43. Wing/body configuration 
details used for the results presented 
in Fig. 44 (ref. 2). 


ARA FULL 

— POTENTIAL CODE y “ ^3*7 o « 4.2° 



(a) y s 0.37. 


Fig. 44. Pressure coefficient compar- 
isons for the configuration of fig- 
ure 43, = 0. 86 (ref. 2) . 
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Fig. 44. Concluded. 
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0«<aNALP^_» 
OF POOR QOAtttV 



Fig. 45. Wing: yuai-1 age details used 
for the results presented in fig- 
ures 46 and 47, EAE wing A + body B8 
(ref. 137). 


Fig. 46. Pressure coefficient com- 
parison for the configuration of fig- 
ure 45, = 0.9, a = 2° (ref. 137). 
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* , 

v>^ J(yp^- 



>0 


Fig. 47. The development of lift with 
iteration number for the configuration 
of figure 45, M„ = 0.9, a = 1° 

(ref. 137). 


o TEST DATA (M - 0.79) 
CALCULATION (M •= 0.80) 



Fig. 48. Comparison of FL022 results 
with experiment for W45j4 supercritical 
wing. Re = 2.4 million (ref. 195). 


O TEST DATA |M - 0.790) 
CALCULATION (M - 0.797) 



Fig. 49. CampariBon of FL027 results 
with experiment for NASA supercritical 
wing. Re = 2.4 million (ref. 195). 


104 













PAGE IS 
W POOR OUAUIY 



Fig, 51. Boeing 747-200 wing/ fuselage 
aonfigui'at ion used ui obtain the 
results of figuPi; 52 (ref. 8), 


ft 


OO EXPERIMENT 



(bJ 


OD EXPERIMENT 




(aJ I'c; 

Fig. 52. Pressure coefficient comparisons for the configuration of figure 51, 
Mg, 0.84, a - 2.8°, with boundary layer (ref. 2). 
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Fig. 5S. Pressure aoeffiaient comparisons for the ONERA M6 wing, M„ 
a = 3.06° (ref. 133). 
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^ig. 54. Comparison of shock sonic Fig. 55. Convergence characteristics 
line positions for the OEERA M6 wing. Of TWING for the ONERA M6 wing, 
yi„ = 0.84, a = 3.06° (ref. 133). =•• 0.84, a = 3.06° (ref. 133). 






’ k/c 


Fig. 56. Comparison of partially and fully converged pressure coefficient 
distributions at ri = 0.20, ONERA M6 wing, = 0.64, a = 3.06° (ref. 133). 


"”oo 

A 7o.8S B.9 experiment 



Fig. 57. Comparison of WING pressure distribution with experiment, Hinson- 
Burdges Ning C (ref. 191). 
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(a) n s 0.3. 
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r x/c 

(b) X] S 0.7. 

Fig. 58. Effect of Mach number and 
angle-of-attaok corrections on the 
Hinson-Burdges Wing C pressure coeffi 
dent distribution (ref. 191). 
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Fig. 59. Compccrison of FL030 pressure 
coefficient distribution with experi- 
ment, three-dimensional viscous cor- 
rections including wake with thickness 
and curvature effects modeled, Hinson- 
Burdges Wing A, = 0.819, a = 1.96° 
(ref. 199). 
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Fig. 60. Lift reduction owing to vis- 
cous effects, R, versus semispan sta- 
tion, advanced transport configuration 
of reference 196 (ref. 199). 


